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Summary 


The  atomization  of  liquid  jets  is  a  fundamental  problem  of  particular  interest  in  determining  the 
performance  and  combustion  stability  of  liquid  rocket  engines.  This  research  effort  is  focused  on  a 
new  analytic  treatment  of  this  free  surface  flow  utilizing  Boundary  Element  Methods  (BEMs).  This 
analytic  approach  is  attractive  for  problems  of  this  nature  due  to  its  inherent  efficiency  (as  com¬ 
pared  to  more  traditional  CFD  methods)  and  its  high  accuracy  in  determining  the  time-dependent 
evolution  of  the  free  surface  for  a  reasonable  computational  expense. 

Ultimately,  the  goal  of  this  project  is  aimed  at  developing  a  complete  simulation  of  an  atom¬ 
ization  process.  Toward  this  end,  we  have  completed  an  axisymmetric,  inviscid  simulation.  This 
step  was  necessary  in  allowing  us  to  develop  a  general  surface  tracking  capability.  In  connection 
with  this  approach,  a  new  capability  has  been  developed  to  model  gas/liquid  flows  with  a  free 
surface  utilizing  the  BEM  approach.  In  addition,  a  new  technique  has  been  developed  to  model  an 
atomization  process  well  beyond  the  initial  droplet  shedding  event  for  an  arbitrary  orifice  geometry 
and  injection  conditions. 

1  Research  Objectives 

The  understanding  of  the  complex  combustion  phenomena  present  in  liquid  rocket  engines  begins 
with  the  fundamental  process  of  fuel  and  oxidizer  jet  atomization.  Since  the  atomization  process  can 
be  greatly  effected  by  acoustic  disturbances1,  it  appears  as  a  primary  focus2  in  studies  involving 
combustion  stability.  For  this  reason,  both  analytic  and  experimental  research  are  necessary  to 
increase  our  understanding  of  the  complex  interactions  between  the  droplet  field  formed  by  the  jet 
and  the  combustion  process. 

The  analytic/numerical  techniques  described  in  this  report  center  on  the  use  of  Boundary  Ele¬ 
ment  Methods  (BEMs)  which  can  be  used  to  accurately  describe  the  time-dependent  evolution  of 
waves  forming  on  the  surfaces  of  the  jet.  The  BEM  approach  is  superior  to  more  traditional  Compu¬ 
tational  Fluid  Dynamic  (CFD)  methods  since  the  BEM  formulation  requires  a  discretization  of  only 
the  external  surface  of  the  jet  rather  than  the  entire  liquid  domain.  This  distinction  is  important 
for  the  highly-distorted  surface  shapes  (and  hence  highly-distorted  CFD  grids)  required  to  properly 
model  these  processes.  In  addition,  a  complex  regridding  would  be  required  of  a  traditional  CFD 
implementation  as  droplets  are  shed  from  the  periphery  of  the  jet.  Models  developed  under  this 
research  grant  have  successfully  demonstrated  this  capability  utilizing  the  BEM  approach. 

With  these  ideas  in  mind,  this  research  effort  seeks  to  develop  a  general  BEM  capable  of 
modeling  the  atomization  process  for  an  arbitrary  fluid  under  arbitrary  injection  conditions  and 
orifice  geometry.  This  capability  would  be  utilized  to  predict  the  primary  atomization  event  and 
permit  interface  to  secondary  atomization  and  collision/ coalescence  models  further  characterizing 
the  spray.  The  following  section  discusses  the  technical  status  of  the  research.  Section  3  describes 
professional  act  ivies,  while  Section  4  gives  the  references  for  this  report. 

2  Status  of  the  Research 

Due  to  the  unusual  starting  time  for  the  contract  (1  July,  1992),  graduate  student  participation 
in  this  research  effort  was  delayed  for  five  months.  As  a  result  of  this  situation,  a  four  month 
contract  extension  was  requested  (and  granted)  as  indicated  in  the  letters  shown  in  Appendix  A. 
This  development  led  to  a  new  schedule  for  the  grant  reporting  requirements  (see  Appendix  A) 
which  changed  the  Final  Report  Due  Date  from  31  August,  1994  to  31  December,  1994.  This 
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document  therefore  serves  to  report  on  the  entire  28  month  grant  period  dating  from  1  July,  1992 
to  31  October,  1994. 

During  this  timeframe,  two  inviscid  computer  codes  have  been  developed  and  the  development 
of  a  viscous  code  is  about  50%  complete. 

The  first  of  the  inviscid  codes  demonstrated  a  capability  of  modeling  an  atomization  process 
for  an  orifice  of  arbitrary  geometry  through  a  large  number  of  droplet  pinchoff  events.  The  model 
represents  the  first  numerical  approach  which  includes  the  impact  of  orifice  geometry  (i.e.  finite  jet 
length)  on  the  atomization  process  within  the  flow  regimes  of  interest.  This  model  is  described  in 
Section  2.1  of  this  report.  The  second  inviscid  model  represents  the  genesis  of  a  new  BEM  capability 
for  general  gas/liquid  flows  involving  a  free  surface.  The  model  has  been  utilized  to  characterize 
atomization  of  a  liquid  jet  in  the  wind-induced  atomization  regime.  This  model  is  described  in 
Section  2.2  of  this  document.  Development  of  a  viscous  capability  are  summarized  in  Section  2.3, 
while  future  efforts  to  be  applied  to  this  research  area  are  discussed  in  Section  2.4. 

2.1  BEM  for  a  Finite-Length  Liquid  Jet 

Historically,  most  analytic  or  numerical  liquid  jet  simulations3*”10  have  ignored  the  presence  of  the 
orifice  and  its  contribution  to  the  atomization  process.  This  assumption  is  quite  valid  in  the  event 
that  the  breakup  of  the  jet  is  far  downstream  of  the  injection  point.  However,  in  the  atomization 
regime  of  interest  to  liquid  rocket  engine  injectors,  jet  breakup  occurs  in  the  immediate  proximity 
of  the  orifice.  For  this  reason,  it  is  necessary  to  include  the  presence  of  the  orifice  (and  its  specific 
design  geometry)  in  any  simulation  within  this  regime. 

In  addition,  to  develop  a  complete  model  of  the  spray  formation  process,  it  is  also  necessary 
to  run  a  simulation  through  numerous  droplet  pinchoff  events.  While  this  capability  has  been 
demonstrated  by  researchers  utilizing  BEMs,  it  has  generally  been  included  to  address  the  dynamics 
of  a  single  droplet  shed  from  a  column  of  fluid9,10.  Here,  a  more  general  capability  is  desired  such 
that  the  dynamics  of  the  main  body  of  fluid  can  be  addressed  after  the  pinchoff  event.  Appendix  B 
contains  a  manuscript  recently  submitted  to  Atomization  and  Sprays  which  describes  the  elements 
of  the  model  created  for  simulation  of  a  finite-length  liquid  jet.  Additional  details  can  be  found  in 
Refs.  11  and  12.  A  brief  summary  of  the  elements  of  this  model  is  provided  below. 

Figure  1  (also  see  Appendix  B,  Fig.  5)  shows  the  computational  domain  utilized  for  a  typical 
simulation.  A  hybrid  grid  is  employed  which  enables  nodes  within  the  orifice  passage  to  remain 
fixed,  while  nodes  on  the  free  surface  are  tracked  in  a  Lagrangian  sense  with  the  local  fluid  velocity. 
Axisymmetric,  inviscid  flow  is  assumed  within  the  liquid,  and  pressure  variations  in  the  gas  (due  to 
the  presence  of  the  liquid)  are  neglected.  These  simplifying  assumptions  permit  us  to  focus  on  the 
development  of  a  scheme  for  tracking  the  complex  evolution  of  a  free  surface  during  atomization 
processes.  A  precise  definition  of  the  free  surface  is  required  for  flows  of  this  nature  in  which  surface 
tension  (and  hence  surface  curvature)  plays  a  substantial  role. 

The  major  features  incorporated  into  this  model  include: 

•  Laplace  Solver  -  This  module  uses  the  BEM  formulation  to  solve  the  differential  equation 
(Laplace’s  equation)  in  the  liquid  domain.  If  <j>  represents  the  velocity  potential  in  the  liquid, 
we  are  focusing  on  the  solution  of  V2<£  —  0  in  this  module. 

•  Time  Integration  -  A  high  accuracy,  fourth-order  time  integration  scheme  (using  the  Runge- 
Kutta  method)  is  employed  to  insure  precise  movement  of  the  surface  nodes.  The  Bernoulli  and 
kinematic  boundary  conditions  specify  the  time-dependent  behavior  of  the  velocity  potential 
and  surface  velocity.  These  equations  are  integrated  using  the  Runge-Kutta  approach.  The 
time  step  for  the  integration  procedure  is  set  dynamically  during  the  integration  in  order  to 
insure  that  nodes  move  only  a  small  fraction  of  the  local  nodal  spacing  on  a  given  time  step. 


2 


Figure  1:  Finite-Length  Liquid  Jet  Computational  Domain. 


•  Free  Surface  Module  -  Fourth-order  centered  differencing  is  used  to  accurately  define  surface 
curvature  at  all  times.  Surface  slope  is  obtained  using  the  method  of  Medina,  et.  al.13. 

•  Surface  Regridding  Module  -  To  avoid  “bunching”  of  points  in  high  curvature  regions,  nodes 
are  redistributed  at  equal  length  intervals  along  the  surface  at  each  time  step.  This  regridding 
procedure  is  implemented  by  fitting  the  surface  coordinates  (and  the  velocity  potential)  with 
cubic  splines.  Points  are  then  redistributed  at  equal  length  intervals  by  interpolating  values 
from  cubic  spline  curvefits.  In  addition,  this  module  includes  the  capability  to  dynamically 
change  the  number  of  nodes  along  the  surface.  Nodes  can  be  added  near  the  orifice  as  more 
fluid  issues  from  the  nozzle.  In  addition,  nodes  can  be  removed  from  the  simulation  when 
droplet  pinchoff  occurs. 

The  Time  Integration,  Free  Surface,  and  Surface  Regridding  Modules  are  developments  which 
will  be  common  to  future  BEM  codes  which  will  include  viscous  effects.  For  this  reason,  it  is  obvious 
that  the  development  of  an  inviscid  model  provides  great  insight  and  capability  to  tackle  the  more 
difficult  viscous  problem.  Numerous  results  of  simulations  performed  with  the  model  are  included 
in  Appendix  B.  Of  particular  interest  is  the  simulation  of  a  dripping  faucet.  This  simulation  is  quite 
appropriate  in  that  atomization  occurs  in  close  proximity  to  the  orifice,  thus  providing  a  means  to 
validate  this  aspect  of  the  model. 

The  geometry  shown  in  Fig.  1  was  used  to  simulate  the  experiments  of  Wu  and  Schelly14. 
These  experiments  revealed  both  chaotic  and  bimodal  droplet  size  formation  within  the  dripping 
flow  regime.  These  authors  attempted  to  correlate  this  behavior  with  chaos  theory.  Results  of  our 
simulation  are  shown  in  Fig.  2  below.  The  results  indicate  that  both  chaotic  and  bimodal  droplet 
size  formation  is  predicted  using  the  deterministic  model.  In  analyzing  results,  we  found  that  this 
behavior  is  attributed  to  interactions  of  disturbances  (caused  by  the  pinchoff  event)  with  the  orifice. 
Therefore,  the  simulation  lends  confidence  to  a  similar  treatment  which  will  be  applied  to  a  viscous 
simulation  in  the  coming  year. 

2.2  Atomization  Processes  in  the  Wind-Induced  Regime 

Atomization  processes  inside  liquid  rocket  engines  can  be  greatly  effected  by  the  influence  of  high 
density  combustion  gases  influencing  the  pressure  distribution  on  liquid  ligaments  formed  during 


3 


2-5. 


1  5h 


1 8  , 5 

5  !»' 


*  i  I 


g  O 


0 

!  8 
0  8 


§  8 


•  • 


0  0  • 


0.4  0.6  0.8  1  \Z  1.4  1.6  1.8 

Weber  Number 


Figure  2:  Droplet  Size  for  Dripping  Flow,  Bo=0.204. 


early  stages  of  atomization.  To  address  the  effects  of  atomization  in  this  “wind-induced”  regime,  a 
BEM  simulation  was  developed15-17  in  accordance  with  the  developments  described  above.  A  de¬ 
tailed  description  of  the  model  and  results  of  simulations  is  attached  in  Appendix  C.  Brief  highlights 
of  the  results  will  be  discussed  below. 

For  this  calculation,  an  infinite  jet  with  a  very  small  perturbation  is  assumed  for  modeling 
purposes.  In  accordance  with  these  assumptions,  the  computational  domain  is  shown  in  Fig.  3 
(adopted  from  Fig.  2  in  Appendix  C).  The  boundary  element  formulation  is  utilized  to  solve  for 
both  liquid  and  gas  flowfields.  The  major  development  required  for  a  model  of  this  type  is  a  stable, 
accurate  treatment  of  the  coupled  set  of  nonlinear  boundary  conditions  present  at  the  gas/liquid 
interface.  In  Appendix  C,  the  new  treatment  developed  for  this  application  is  described  in  detail. 
Results  indicate  that  the  methodology  employed  is  quite  general,  and  can  easily  be  applied  to  a 
variety  of  other  gas/liquid  flow  problems. 

The  model  has  been  shown  to  reproduce  both  linear  and  nonlinear  results  of  previous  researchers. 
Moreover,  the  model  predicts  the  presence  (and  size)  of  satellite  droplets  formed  within  the  first 
wind-induced  regime.  Figure  4  depicts  a  comparison  of  predicted  droplet  sizes  with  measurements 
of  researchers  indicating  a  good  predictive  capability.  The  model  has  also  been  utilized  to  predict 
the  transition  between  first  and  second  wind-induced  regimes.  This  transition,  which  can  only 
be  evaluated  with  a  nonlinear  model,  corresponds  to  the  point  where  droplets  are  shed  at  the  jet 
periphery,  rather  than  at  the  centerline  of  the  jet.  Figure  5  presents  results  of  a  simulation  near 
this  transition  point.  Results  indicate  that  the  transition  occurs  at  a  gas-based  Weber  number  of 
approximately  2.5. 

In  addition,  the  code  has  been  utilized  to  model  the  development  of  ripples  on  the  surface  of  a 
high-speed  liquid  jet  operating  in  the  second  wind-induced  regime.  Here  we  should  note  that  this 
regime  is  consistent  with  injection  conditions  in  some  liquid  rocket  engine  injectors.  A  simulation 
of  the  process  is  shown  in  Figure  6,  The  inviscid  gas-phase  formulation  employed  leads  to  a 
“spiked”  appearance.  Inclusion  of  gas-phase  viscosity  is  required  to  improve  the  simulation  at  large 
distortions  (where  boundary  layer  separation  are  presumed  to  occur).  Nevertheless,  this  model 
represents  the  first  nonlinear  simulation  of  gas/liquid  flows  within  this  high  speed  flow  regime. 
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Figure  3:  Schematic  of  Computational  Domain  Denoting  Boundary  Conditions 
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Figure  4:  Comparison  of  the  BEM  to  the  Drop  Sizes  of  Rutland  and  Jameson8  (o)  and  Lafrance 
(x)  for  Main  and  Satellite  Drops,  e  =  0.00129 
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Figure  5:  Nonlinear  Jet  Evolution  in  the  Transition  Region 


Figure  6:  Nonlinear  Jet  Evolution  in  the  Second  Wind-Induced  Regime 
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2*3  Current  Status  of  Viscous  Code  Development 

Since  this  research  project  has  been  granted  a  continuation  (AFOSR  Contract  F49620-94-1-0151), 
it  is  appropriate  to  discuss  the  future  directions  of  this  work  within  the  context  of  this  final  report. 
As  noted  in  the  previous  discussion,  the  ultimate  goal  of  this  work  is  to  establish  a  simulation  of 
the  viscous  atomization  process.  To  this  end,  we  began  the  work  based  on  the  simpler  inviscid 
formulation  in  order  to  gain  experience  in  boundary  element  techniques  and  to  validate  free  surface 
modules  required  for  viscous  simulations. 

Prior  to  the  development  of  a  viscous  BEM  capability,  an  extensive  literature  review  was  con¬ 
ducted.  While  several  authors  have  developed  viscous  simulations  using  BEM  techniques,  many 
of  these  methods  use  alternate  approaches  in  viscous  zones18,19,  or  are  restricted  to  Stokes  flow 
regimes  at  very  low  Reynolds  numbers20"22.  Recent  developments23-25  indicate  that  a  full  viscous 
BEM  solution  is  now  at  hand.  In  particular,  the  methodology  developed  by  Jin  and  Brown25,  is 
very  attractive  for  free  surface  problems  in  that  it  employs  a  Lagrangian  point  of  view  which  is 
natural  to  use  for  tracking  a  free  surface. 

Jin  and  Brown  described  a  method  for  treating  the  Lagrangian  time  derivative  using  the 
traction-velocity  boundary  element  formulation.  This  formulation  has  been  applied  successfully  to 
steady  free  surface  flows  by  various  researchers26,27.  It  is  important  to  note  that  Jin  and  Brown’s 
technique  for  unsteady  problems  has  never  been  implemented  in  an  actual  code. 

For  the  unsteady  problem,  Jin  and  Brown  propose  a  treatment  utilizing  the  Dual  Reciprocity 
Method ,  (DRM)  to  aid  in  mapping  domain  integrals  (which  appear  due  to  presence  of  unsteady 
terms)  to  the  boundary.  The  DRM  represents  a  new  technique  developed  within  the  past  decade 
by  the  BEM  research  group23  at  Wessex  Institute  of  Technology  in  Southampton,  United  Kingdom. 
Dr.  Carlos  Brebbia,  the  world’s  leading  authority  on  BEM  development,  heads  this  research  group. 
It  is  only  within  the  past  two  years  that  a  book  on  DRM  implementation  has  become  available, 
effectively  maps  domain  integrals  to  the  boundary  by  utilizing  shape  functions  which  are  forms  of  the 
governing  equations  themselves.  Only  a  few  interior  nodes  must  be  utilized  in  such  a  formulation; 
an  attractive  feature  for  boundary-based  solution  techniques. 

To  develop  a  free  surface  code  utilizing  the  methodology  described  above  is  quite  daunting  since 
a  full  unsteady  capability  using  the  DRM  is  required.  For  this  reason,  an  incremental  approach  was 
adopted  whereby  a  series  of  models  with  increasing  complexity  are  developed.  While  the  ultimate 
goal  is  to  create  an  axisymmetric,  viscous,  unsteady  code,  the  following  approach  has  been  adopted: 

1.  Develop  a  steady,  two-dimensional  viscous  code.  Validate  this  model  by  performing  a 
simulation  of  general  Couette  flow  (flow  in  a  channel)  using  various  computational  do¬ 
mains.  This  code  provides  insight  into  the  tract  ion- velocity  formulation  and  the  required 
treatment  of  singularities  and  multi-valued  tractions  at  domain  corners. 

2.  Test  case  simulation  using  the  DRM  on  a  problem  with  a  known  analytic  solution28. 
This  validation  aids  in  gaining  experience  in  applying  the  DRM  to  complex  differential 
equations. 

3.  Extend  the  two-dimensional  steady  model  to  handle  unsteady  problems  via  implementa¬ 
tion  of  the  DRM  on  the  system  of  viscous  equations.  Validate  this  code  by  performing  a 
simulation  of  Stokes’s  first  problem  (impulsively-started  infinite  flat  plate).  This  develop¬ 
ment  will  permit  the  determination  of  the  number  of  internal  nodes  required  for  accuracy 
with  the  DRM. 

4.  Apply  the  Lagrangian  treatment  of  boundary  nodes  per  Jin  and  Brown  to  the  two- 
dimensional,  unsteady  code.  The  addition  of  a  free  surface  module  at  this  point  will  give 
the  capability  to  investigate  the  stability  of  liquid  sheets. 
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5.  Development  of  an  axisymmetric  model  using  the  methodology  created  for  the  two- 
dimensional  code.  A  steady  model  will  verify  the  new  integrations  in  the  traction-velocity 
formulation,  a  second  model  will  verify  the  axisymmetric  formulation  of  the  DRM,  fol¬ 
lowed  by  a  full  unsteady  model  in  Lagrangian  form. 

At  the  present  time,  Tasks  #1  and  2  have  been  completed,  and  Task  #3  is  well  underway.  The 
two-dimensional  formulation  was  implemented  first  due  to  the  relative  simplicity  of  the  fundamental 
solutions  (Green's  functions)  for  this  case.  For  the  axisymmetric  case,  the  Green's  functions  are 
much  more  complex.  Items  #4  and  5  will  be  accomplished  under  a  follow-on  contract  which  will 
be  concluded  in  June,  1996. 

2.4  Future  Efforts  -  Viscous  Flows 

Upon  completion  of  the  tasks  outlined  in  the  previous  list,  we  shall  have  the  capability  of  modeling 
both  two-dimensional  and  axisymmetric  flows  with  free  surfaces.  These  new  models  will  aid  in 
enhancing  the  understanding  of  complex,  nonlinear  atomization  processes  germaine  to  both  rocket 
and  airbreathing  propulsion  systems.  In  particular,  two-dimensional  simulations  of  liquid  sheets 
should  be  underway  within  the  next  two  months.  Another  interesting  two-dimensional  problem 
involves  the  interaction  of  a  liquid  sheet  with  a  transverse  gas  flow.  This  flowfield  is  important  to 
atomization  processes  in  current  state-of-the-art  gas  turbine  engine  fuel  injectors.  To  perform  this 
calculation,  gas  flow  interactions  (such  as  those  discussed  in  Section  2.2)  will  need  to  be  incorporated 
into  the  model. 

Upon  the  completion  of  the  axisymmetric  code,  a  variety  of  droplet  and  liquid  jet  simulations 
can  be  conducted.  Droplet-related  problems  such  as  secondary  atomization  and  droplet  colli¬ 
sion/coalescence  will  be  considered.  Blanching  problems  associated  with  liquid  oxygen  droplet 
impact  on  chamber  walls  can  also  be  investigated.  In  addition,  a  complete  simulation  of  the  liquid 
jet  (including  internal  orifice  flow)  will  be  at  hand  with  the  completion  of  the  axisymmetric  model. 
This  simulation  will  be  useful  in  gaining  knowledge  of  primary  atomization  processes  which  are 
difficult  (if  not  impossible)  to  observe  experimentally  due  to  the  obscuration  of  this  region  by  the 
remainder  of  the  spray. 

3  Professional  Activities 

Currently  there  is  one  Ph.D.  student,  Mr.  J.  H.  Hilbing,  working  on  this  project.  Three  masters 
students;  Mr.  I.  Murray,  Mr.  M.  Moses,  and  Mr.  M.  Rutz  will  be  performing  related  research 
through  funding  received  in  a  recently-awarded  Augmentation  Award  for  Science  and  Engineering 
Research  Training  (AASERT).  As  a  result  of  this  sponsorship,  Mr.  C.  Spangler  recently  received 
his  Masters  Degree: 


Spangler,  C.A.,  “Nonlinear  Modeling  of  Jet  Breakup  in  the  Wind-Induced  Regime,”  M.S.  Thesis, 
Purdue  University,  1994. 

At  the  present  time,  two  publications  (attached  in  Appendices  B  and  C)  have  been  submitted 
to  referred  journals  as  a  result  of  this  research: 


Hilbing,  J.  H.,  and  Heister,  S.  D.,  “A  Boundary  Element  Method  for  Atomization  of  a  Finite 
Liquid  Jet”,  Submitted  to  Atomization  and  Sprays ,  September,  1994. 
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Spangler,  C.  A.,  Hilbing,  J.  H.,  and  Heister,  S.  D.,  “Nonlinear  Modeling  of  Jet  Atomization  in  the 
Wind-Induced  Regime”,  Submitted  to  Physics  of  Fluids,  August,  1994. 

In  addition,  the  authors  of  these  papers  attended  this  year’s  conference  sponsored  by  the  In¬ 
stitute  for  Liquid  Atomization  and  Spray  Systems  (ILASS),  which  was  held  in  Seattle,  WA  during 
June,  1994.  Both  works  were  presented  at  this  forum;  in  addition,  Professor  Heister  served  as  a 
member  of  a  Panel  Discussion  on  Spray  Modeling.  These  items  are  summarized  below: 


Hilbing,  J.H.,  and  Heister,,  S.D.,  UA  Boundary  Element  Method  for  Liquid  Jet  Atomization 
Processes”,  ILASS-94  Conference  Proceedings,  June  1994. 

Spangler,  C.  A.,  and  Heister,  S.  D.,  “Nonlinear  Modeling  of  Jet  Atomization  in  the  Wind-Induced 
Regime”,  ILASS-94  Conference  Proceedings,  June,  1994. 

Heister,  S.  D.,  “Using  Boundary  Element  Methods  to  Model  Atomization  Processes”,  Invited 
Lecture/  Panel  Member,  ILASS-94  Tutorial  Workshop  on  Spray  Modeling,  June,  1994. 

Other  related  work,  sponsored  by  Cummins  Engine  Company,  has  been  developed  during  the 
past  year: 


Chen,  Y  and  Heister,  S.  D.  A  numerical  treatment  for  attached  cavitation,  to  appear  in  J.  Fluids 
Engineering  1994. 

Chen,  Y.  and  Heister,  S.  D.  Two-Phase  Modeling  for  Cavitating  Flows.  Cavitation  and 
Gas-Liquid  Flows  in  Fluid  Machinery  and  Devices.  FED-Vol.  190  299-307,  1994. 

Chen,  Y.,  and  Heister,  S.  D.,  “Modeling  Hydrodynamic  Non-Equilibrium  in  Bubbly  and 
Cavitating  Flows”,  Submitted  to  Journal  of  Fluid  Mechanics ,  1994. 

We  are  currently  hoping  to  use  these  models  to  address  cavitating  orifice  experiments  recently 
conducted  at  Phillips  Lab  under  the  supervision  of  Dr.  Doug  Talley. 

Finally,  another  manuscript  resulting  from  funding  of  this  research  is  currently  in  preparation: 


Heister,  S.  D.,  “Modeling  Gas/Liquid  Flows  using  Boudary  Element  Methods”,  to  be  submitted, 
International  Journal  of  Numerical  Methods  in  Fluids ,  1994. 

3.1  Technology  Transfer/Coupling  Activities 

This  research  is  rapidly  approaching  the  point  at  which  significant  technology  transfer  and  cou¬ 
pling  activities  will  occur.  While  the  goal  of  this  research  lies  in  improving  the  understanding 
of  atomization  processes  in  connection  with  combustion  instability  in  liquid  rocket  engines,  the 
models  being  developed  are  quite  general  in  nature  and  can  be  applied  to  many  other  problems. 
Examples  of  commercial  applications  of  interest  include:  paint  sprays,  polymer/wax  extrusion  and 
manufacturing  processes,  agricultural  sprays,  diesel  engine  injectors,  inkjet  printers,  and  fire  hoses 
to  name  just  a  few. 

There  was  a  significant  interest  in  these  modeling  efforts  at  the  ILASS  meeting  attended  during 
June.  Not  only  was  the  principle  investigator  invited  to  serve  as  a  panel  member  for  a  tutorial 
workshop  (see  above),  but  numerous  industrial  officials  expressed  a  strong  interest  in  this  work.  As 
a  result,  we  have  developed  coupling  activities  with  Mr.  C.  Lipp  of  Dow  Chemical  Company  and 
Mr.  G.  P.  Anath  of  Johnson  Wax  Corporation.  Both  of  these  individuals  are  eager  to  apply  our 
future  viscous  modeling  capabilities  to  industrial  manufacturing  processes.  Officials  at  Cummins 
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Engine  Company  (Dr.  W.  Eckerle  and  Dr.  S.  Menon)  are  also  interested  in  developments  of  this 
modeling  capability  in  connection  with  injectors  for  diesel  engines.  Individuals  involved  with  this 
research  are  excited  about  the  industrial  applications  of  this  capability  and  look  forward  to  possible 
interactions  with  the  industrial  officials  mentioned  above. 
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Nonlinear  modeling  of  jet  atomization  in  the  wind-induced  regime 

Christopher  A.  Spangler,  James  H.  Hilbing,  and  Stephen  D.  Heister 

Department  of  Aeronautics  and  Astronautics,  Purdue  University,  West  Lafayette,  Indiana  47906 
(Received 


Abstract 

A  Boundary  Element  Method  (BEM)  has  been  developed  to  solve  for  the  nonlinear  evolution  of  a  liquid 
jet  acting  under  the  influence  of  both  surface  tension  and  the  aerodynamic  interactions  with  the  surrounding 
atmosphere.  For  longer  waves,  aerodynamic  effects  are  shown  to  cause  a  “swelling”  of  the  liquid  surface  in 
the  trough  region.  The  model  predicts  the  presence  of  satellite  drops  in  the  first  wind-induced  regime,  and 
predicts  the  evolution  of  a  “spiked”  surface  at  the  periphery  of  the  jet  for  conditions  consistent  with  the 
second  wind-induced  regime.  The  effects  of  the  disturbance  wave  number,  the  liquid  Weber  number,  and  the 
density  ratio  between  the  liquid  jet  and  the  surrounding  gas  on  the  breakup  of  the  jet  have  been  examined. 
Transition  points  between  various  flow  regimes  have  also  been  identified. 
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I.  Introduction 


The  atomization  of  a  liquid  jet  in  a  gaseous  environment  is  a  fundamental  problem  in  two-phase  flows. 
The  numerous  applications  of  this  flow  include  devices  used  for  liquid  fuel  injection  systems,  spray  painting, 
and  ink  jet  printing.  Reitz  and  Bracco1  developed  a  classification  of  liquid  jet  atomization  regimes  as 
highlighted  in  Fig.  1.  The  Rayleigh  and  first  wind-induced  regimes  are  characterized  by  the  formation  of 
drops  with  diameters  larger  than  that  of  the  orifice.  The  breakup  in  the  Rayleigh  regime  is  governed  solely  by 
capillary  instabilities,  while  in  the  first  wind-induced  regime,  the  capillary  instabilities  are  enhanced  by  the 
aerodynamic  interactions  with  the  surrounding  gas.  For  the  second  wind-induced  and  atomization  regimes, 
the  breakup  is  characterized  by  the  formation  of  drops  significantly  smaller  than  the  jet  diameter.  In  this 
case,  the  wave  formation  is  driven  by  aerodynamic  forces  since  surface  tension  serves  as  a  stabilizing  force 
for  the  shorter  waves  formed  in  this  regime.  In  the  atomization  regime,  the  drops  are  formed  very  near  the 
orifice  exit  plane  and  viscous  effects  are  known  to  be  of  great  importance. 

The  theoretical  basis  for  the  capillary  instability  was  first  investigated  by  Rayleigh,2  who  conducted  a 
linear  analysis  neglecting  the  effect  of  gas  pressure  variations  on  jet  distortion.  Weber3  developed  a  linear 
theory  similar  to  Rayleigh’s,  but  he  included  the  effects  of  both  the  liquid  viscosity  and  the  pressure  of  the 
surrounding  gas  on  the  jet  behavior.  The  effect  of  variations  in  gas  velocity  was  later  incorporated  into  a 
linear  analysis  by  Sterling  and  Sleicher.4  By  considering  higher-order  terms  in  an  expansion  of  the  surface 
deformation  or  jet  velocity,  more  recent  efforts5”7  have  been  able  to  improve  on  predicted  surface  distortions 
due  to  the  capillary  instabilities  in  the  absence  of  gas-phase  pressure  variations.  These  works  predicted  that 
a  single  wave  would  divide  into  two  droplets  (called  main  and  satellite  drops)  as  opposed  to  the  single  drop 
predicted  by  the  linear  theory. 

Experimental  measurements  of  flows  consistent  with  the  Rayleigh  regime7-10  confirmed  the  presence  of 
satellite  drops  formed  during  the  nonlinear  portion  of  the  atomization  process.  These  experiments  were  per¬ 
formed  by  using  a  forced  oscillation  frequency  thus  permitting  investigation  of  the  growth  of  instabilities  over 
a  range  of  wavelengths.  More  recently,  experimental  efforts  have  also  focused  on  the  control  of  the  satellite 
size  through  the  use  of  higher  harmonics  of  the  forcing  frequency11,12  as  well  as  amplitude  modulation13 
techniques. 

Unfortunately,  very  few  experiments  have  been  reported  at  conditions  consistent  with  the  first  wind- 
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induced  regime.  While  some  efforts  have  focused  on  jet  breakup  length,4’9  very  few  droplet  size  measurements 
have  been  made  at  moderate  jet  velocities  consistent  with  the  first  wind-induced  regime.  The  influence  of 
viscosity  on  the  breakup  of  a  liquid  jet  was  examined  by  Goedde  and  Yuen.10  They  found  that  viscosity 
served  as  a  stabilizing  influence  on  the  growth  rate.  They  also  noted  that  the  pressure  field  around  the  liquid 
jet  has  an  influence  on  the  formation  of  droplets  although  they  did  not  quantify  this  effect.  Rutland  and 
Jameson14  examined  the  occurrences  of  secondary  swellings  between  the  primary  crests  of  a  liquid  jet. 

Most  recently,  Mansour  and  Lundgren15  developed  a  full  nonlinear  simulation  of  the  Rayleigh  breakup 
process  utilizing  a  boundary  element  method  (BEM).  The  use  of  a  BEM  is  attractive  in  flows  of  this  nature 
because  of  the  high  accuracy  these  techniques  provide  in  assessing  curvature  of  highly  deformed  surfaces. 
This  feature  is  crucial  in  flows  in  which  capillary  forces  play  a  substantial  role.  This  work  presents  an 
extension  of  the  model  due  to  Mansour  and  Lundgren  for  the  case  where  gas-phase  pressure  variations  are 
not  negligible.  The  model,  based  on  a  BEM  approach,  is  discussed  in  the  following  section.  Comparisons  are 
made  for  both  linear  theory  and  measured  drop  sizes,  and  detailed  surface  evolutions  are  presented  as  well. 

II.  Model  Development 

In  developing  the  model,  we  assume  that  the  flow  is  inviscid  and  incompressible,  and  that  no  gravity  or 
other  body  forces  are  present.  Under  these  assumptions,  the  unsteady  liquid  and  gas  behavior  is  described 
by  Laplace’s  equation  V2<d  =  0  where  <j>  is  the  velocity  potential  such  that  its  gradient  is  simply  the  velocity, 
V  •  <j>  —  v.  We  also  assume  that  finite  disturbances  of  a  periodic  nature  develop  well  downstream  of  the  jet 
orifice  such  that  the  evolution  of  the  jet  is  accurately  represented  as  an  infinite  wavetrain.  This  assumption 
removes  the  influence  of  the  orifice  on  the  evolution  of  the  jet  surface. 

For  this  analysis,  the  characteristic  dimensions  utilized  are  the  unperturbed  jet  radius,  the  mean  jet 
velocity,  and  the  liquid  density.  Under  this  nondimensionalization,  the  governing  equations  indicate  that 
the  Weber  number  and  gas/liquid  density  ratio  are  the  appropriate  dimensionless  variables  characterizing 
this  flow.  We  define  the  Weber  number  as:  We  =  piU2a/<r ,  where  p  is  the  liquid  density,  U  is  the  relative 
velocity  between  the  liquid  and  gas,  a  is  the  orifice  radius,  and  <r  is  the  surface  tension.  Also,  the  density 
ratio  is:  c  =  pg/pt  where  pg  is  the  gas  density.  The  nondimensional  form  of  the  Bernoulli  equation  relates 
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velocity  potentials  to  the  local  pressure  at  the  interface.  In  the  liquid  region,  this  condition  is  expressed: 


(1) 


while  in  the  gas  region  it  is: 

^  +  5(W,)2  +  P,=°  (2) 

where  k  is  the  surface  curvature  and  Pg  is  the  gas  pressure  on  the  surface.  These  two  equations  are  boundary 
conditions  on  the  interface  and  are  coupled  through  the  local  gas  pressure.  This  coupling  must  be  addressed 
in  the  development  of  a  consistent,  stable  integration  procedure. 

Figure  2  depicts  the  computational  domain,  coordinate  system,  and  boundary  conditions  utilized  in  the 
calculations.  In  this  figure,  4>g  represents  the  velocity  potential  in  the  gas  phase,  q  represents  the  velocity 
normal  to  the  local  surface,  and  the  radial  and  axial  coordinates  are  denoted  r,  z,  respectively.  Nodes  are 
placed  along  the  interface,  along  vertical  planes  at  both  ends  of  the  wave,  and  along  a  horizontal  surface,  in 
the  gas,  far  from  the  interface.  The  liquid  nodes  are  denoted  “o”,  while  the  gas  nodes  are  denoted  Mx”.  The 
initial  wave  shape  of  the  surface  is  assumed  to  be: 


r  =  1  -f  rjcos(kz)  (3) 

where  tj  is  the  initial  surface  deflection  and  k  is  the  wave  number  of  the  forcing  disturbance  on  the  jet.  The 
solution  procedure  employed  requires  an  initial  prediction  for  both  the  pressure,  Pg)  on  the  wave  surface  and 
<f)g  along  the  surface,  which  are  determined  from  the  potential  flow  solution  for  flow  over  a  wavy  cylinder16: 

Pg  = -erik  cos(kz)^i£ll  (4) 

4>g  =  z~V  sin(fcz)^0^  (5) 

where  Ko  and  K\  are  modified  Bessel  functions  of  the  second  kind.  The  ends  of  the  gas  domain  are  gridded 
using  an  exponential  stretching  technique  to  reduce  the  total  number  of  nodes  on  these  surfaces,  while 
retaining  the  higher  accuracy  of  the  finer  mesh  near  the  interface.  The  uppermost  boundary  in  the  gas  phase 
is  place  ten  radii  from  the  centerline.  Numerical  experiments  indicate  that  this  distance  is  adequate  to  apply 
the  q  =  0  boundary  condition  without  adverse  effects  on  the  behavior  of  the  interface.  Nodes  are  initially 
placed  with  equal  spacing  along  all  surfaces  associated  with  the  liquid  domain. 
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The  Boundary  Element  Method  is  used  to  develop  a  computational  model  to  predict  the  nonlinear 
evolution  of  the  jet  surface  in  time17.  The  formulation  of  the  BEM  starts  with  the  integral  representation 
of  Laplace’s  equation.  Following  Liggett  and  Liu18,  this  is: 


^)  +  jf  (6) 

where  rj*  denotes  the  vector  from  the  origin  to  point  “F ,  <^(rj)  is  the  value  of  the  potential  at  this  point, 
T  denotes  the  boundary  of  the  domain,  and  G  is  the  free  space  Green’s  function.  For  the  BEM,  either  <f> 
or  q  =  d<f>/dh  is  specified  at  each  “node”  on  the  boundary.  Here  h  is  the  outward  normal  to  the  boundary 
such  that  q  is  the  velocity  normal  to  the  boundary.  The  quantity  a  in  Eq.  6  results  from  singularities 
introduced  as  the  integration  passes  over  the  boundary  point,  rj.  For  the  axisymmetric  formulation,  a  is 
related  to  the  angle  between  successive  nodes  on  the  boundary.  The  free  space  Green’s  function  solution  to 
the  axisymmetric  Laplacian  is18: 


_  4  rK(p) 


(7) 


with 


a  =  (r  +  ri)2  +  (z  -  z{)2 


(r-r,)2  +  (*-*02 
(r  +  rt)2  +  (*-*,)2 


(8) 


where  K(p)  is  the  complete  elliptic  integral  of  the  first  kind. 

Physically,  the  free  space  Green’s  function  represents  the  influence  of  a  source  at  location  r  of  strength 
27 r  on  the  “base  point”,  t\.  While  the  quantity  dG/dh  represents  the  influence  of  a  doublet  at  location  r 
on  the  base  point.  By  properly  weighting  all  these  source  and  doublet  strengths,  Laplace’s  equation  can 
be  solved  for  arbitrary  boundary  conditions.  Derivation  of  dG/dn  for  this  axisymmetric  problem  is  quite 
involved.  Gipson19  gives: 


dG  _  -1 
dn  2  yja 


nrK(p)  +  (r  _  r.y  +  \z  _  z.)2  (((r  ~  r')2  -  (*  ~  *i )2)  +  2rnz(z  -  z<)] 


(9) 


(r-r,)2  +  (z 

with  the  quantities  nr  and  hz  being  the  components  of  the  outward  normal  unit  vector  in  the  radial  and 
axial  directions,  respectively,  while  E(p)  is  the  complete  elliptic  integral  of  the  second  kind. 

Performing  the  desired  integration  around  the  boundary  of  the  domain  results  in  an  algebraic  system  of 
equations.  In  order  to  perform  this  integration,  a  linear  element  formulation  is  employed  in  which  <j) ,  <j>9 ,  and 
q  are  assumed  to  have  a  linear  variation  between  successive  nodes.  With  this  assumption,  the  values  for  the 
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potential  of  the  element  bounded  by  nodes  j  and  j  +  1  is: 


,  _  ,  ri+i  -  r  r-r, 

v  -  V)  7T— - —  +  <f>j+ 1  i 


(10) 


ri+i-r,-  7  ri+i  -r;- 

with  analogous  expressions  for  <j>3  and  q.  Using  linear  elements,  Eq.  6  is  divided  into  two  separate  integrals: 


L 


rj  —  1 


qGdT  =  (' )  +  qj+ih(i,j)} 

j=i 


and 


where 


K 


n-  1 


*7^  =  E  fo  *»(»'.  j) + tj+iU(ij)i 

i 


_  /r’+‘  ri+1  -  r  /l  j+l  rj+1  -  r  ac 

h-Jr,  iWj""  ,3-/r,  MS' 

_  /~r,+i  r~r>  dGdr 


(11) 

(12) 

(13) 

(14) 


JrL  j+i  T  — r  rl>+i  r_r 

f  r  /4=  /  — — — __e_  , 

r>  r;+i  ~  Vr,  ri+i  -  Ij  <9n 

The  integrals  Ji  through  J4  are  calculated  using  four  point  Gaussian  Quadrature,  and  are  used  to  create 
D  and  5  matrices  which  contain  doublet  and  source  contributions,  respectively.  Since  G  and  dG/dn  are 
singular  when  the  integration  passes  over  the  base  point,  special  treatment  of  the  integrals  in  Eqs.  13  and 
14  is  used  for  these  cases17.  From  Eqs.  11  and  12,  we  can  write: 


Sij  —  Iii  +  hi'j-x 


Di,j  =  hiti  +  hij.t 


(15) 


so  that  the  integration  of  Eq.  6  results  in  a  system  of  algebraic  equations  of  the  form: 


[Dm  =  [S]{q} 


(16) 


with  the  values  of  a  incorporated  into  diagonal  elements  of  the  D  matrix.  To  solve  this  system  of  equations, 
all  the  unknown  terms  are  “pivoted”  by  interchanging  the  columns  of  the  appropriate  matrices,  giving  the 
standard  form  for  a  linear  system  of  equations.  The  system  is  fully  populated  and  is  solved  by  an  LU 
decomposition  method. 

The  boundary  conditions  on  the  interface  between  the  liquid  and  the  gas  are  the  Bernoulli  equations  in 
the  two  domains.  Nodes  on  the  interface  are  assumed  to  travel  with  the  local  liquid  surface  velocity,  so  a 
transformation  from  the  Eulerian  to  Lagrangian  reference  frame  is  required.  After  this  transformation,  Eqs. 
1  and  2  become: 


D<t> 

Dt 


=  \  (W)2  -  p9 


K 

We 


(17) 
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(18) 


^=!-|(Wf)2-e^  +  v*-v*< 

In  these  expressions^  the  notation  D  /Dt  denotes  changes  in  time  for  nodes  moving  with  the  liquid  interfacial 
velocity. 

The  procedure  to  solve  these  two  equations  is  as  follows:  At  the  beginning  of  each  time  step,  an  initial 
value  of  <f>  is  known.  Using  this  value  of  <f>  as  the  boundary  condition  on  the  interface,  solution  of  Laplace’s 
equation  in  the  liquid  domain  provides  values  of  q  on  the  interface.  The  interfacial  velocity  is  then  used  as 
a  boundary  condition  for  the  gas  phase  solution  of  Laplace’s  equation.  Using  this  value  of  g,  the  <j>g  value 
on  the  gas  side  of  the  interface  is  obtained.  The  new  value  of  the  gas  velocity  potential  in  Eq.  18  allows  a 
new  gas  pressure  to  be  calculated  by  approximating  D<f>g/Dt  using  a  first  order  backward  difference  in  time. 
The  new  gas  pressure  is  then  put  into  Eq.  17  to  determine  the  change  in  <f>  with  time. 

Nodes  on  the  surface  are  moved  by  integration  of  the  kinematic  boundary  conditions: 


Dr  _  d<j> 
Dt  dr 


D z  _  d<)> 
Dt  dz 


(19) 


d<t>  d<f>  . 

-  =  -Sinf3  +  qcosf3 
d<t>  d<t> 

-^-COsP-qsmP 


Each  time  <j>  is  updated,  Eqs.  19  are  integrated  to  determine  the  motion  of  nodes  on  the  interface  using  a 
fourth-order  Runge  Kutta  method.  The  velocities  d<j>/dr ,  d<j>/dz  are  obtained  from  a  coordinate  transfor¬ 
mation: 

drft  flrfi 

(20) 
(21) 

Here,  s  is  a  natural  coordinate  aligned  with  the  surface.  The  velocity -tangential  to  the  surface,  d<f>/ds}  is 
obtained  from  a  fourth-order  centered  difference  approximation.  The  quantity  (3  represents  the  local  angle  of 
the  surface  with  respect  to  the  horizontal  direction.  This  angle  is  determined  by  fitting  a  parabola  through 
the  point  in  question  and  its  nearest  neighbors.  Surface  curvatures  are  obtained  by  utilizing  a  fourth-order, 
centered  difference  formula. 

Since  the  nodes  on  the  interface  are  allowed  to  move  with  their  local  velocity,  over  time  they  tend  to 
group  themselves  in  the  center  of  the  jet.  This  grouping  of  the  nodes  is  undesirable  as  it  leaves  the  ends  of 
the  interface  poorly  defined.  To  alleviate  this  problem,  the  surface  mesh  is  regridded  using  a  series  of  cubic 
splines  (for  r,  z,  <j>^  and  4>g)  at  each  time  step  to  keep  the  spacing  between  the  nodes  constant  across  the 
surface.  The  use  of  the  Runge-Kutta  integration  scheme  is  well  suited  to  this  remeshing  scheme,  since  it 
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does  not  require  information  on  node  positions  at  previous  time  levels  to  predict  the  subsequent  motion  of 
the  surface. 

A  special  treatment  is  required  for  nodes  on  each  end  of  the  interface  because  the  outward  normal  vector 
is  multivalued  at  these  corners  in  the  domain.  In  the  liquid  domain,  we  have  q  =  0  in  the  axial  direction 
due  to  the  periodic  boundary  condition,  so  that  the  node  is  tracked  vertically.  In  the  gas  domain,  the  axial 
velocity  is  lagged  one  time  step  and  the  vertical  velocity  is  obtained  from  the  liquid  domain  solution.  Using 
these  two  velocity  components  at  the  corner  permits  the  solution  of  <f>g  at  this  location. 

As  the  interface  deforms  over  time,  eventually  it  will  reach  the  centerline  of  the  domain.  Physically,  at 
this  time,  the  jet  pinches  into  a  series  of  distinct  droplets.  Numerically,  to  prevent  a  node  from  obtaining 
a  negative  radius,  the  jet  is  considered  “pinched”  when  any  node  on  the  surface  comes  within  a  distance  of 
1%  of  the  initial  jet  radius  of  the  centerline.  At  this  time,  the  jet  is  broken  into  discrete  drops,  and  their 
volume  and  radii  are  calculated.  The  radii  of  these  drops  has  been  found  to  be  insensitive  to  this  assumed 
pinch  criteria. 

III.  Model  Validation 

The  model  has  been  validated  through  comparisons  with  previous  experimental  and  computational  stud¬ 
ies.  By  assuming  a  very  small  initial  wave  height,  comparisons  can  be  made  with  the  linear  theory  developed 
by  Sterling  and  Sleicher.4  Figure  3  depicts  the  growth  rate  comparison  as  a  function  of  disturbance  wave 
number  for  different  Weber  numbers  at  a  density  ratio  of  0.001,  approximately  that  of  ambient  air  surround¬ 
ing  a  water  jet.  Errors  of  less  than  1%  are  achieved  using  grids  involving  as  few  as  59  nodes  along  the 
interface.  The  different  Weber  numbers  examined  correspond  to  breakup  in  the  Rayleigh  regime  (We  =  2), 
the  first  wind-induced  regime  (We=  1500),  and  the  second  wind-induced  regime  (We  =  15,000). 

The  initial  surface  deflection  (17)  was  assumed  to  be  0.001  for  these  calculations.  It  is  interesting  to 
note  that  substantial  differences  between  nonlinear  calculations  and  linear  theory  occur  even  for  r)  values  as 
small  as  0.01.  This  result  indicates  that  nonlinear  behavior  can  be  affected  by  relatively  small  perturbations 
induced  by  some  mechanical  device  with  a  minimal  energy  input.  Various  researchers  have  capitalized  on  this 
behavior  in  controlling  droplet  sizes  experimentally  by  using  frequency  modulation11’12  or  through  amplitude 
and  frequency  modulation.13 

Nonlinear  attributes  of  the  calculation  were  validated  by  comparison  with  similar  predictions  of  Mansour 
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and  Lundgren15,  as  well  as  experimental  measurements  of  Rutland  and  Jameson.14  The  comparison  to  the 
nonlinear  surface  shape  calculated  by  Mansour  and  Lundgren15  is  not  presented  here  in  the  interest  of  brevity. 
When  the  current  model  is  run  for  the  same  conditions  (We  =  2,  e  =  0.0,  ib  =  0.3,  0.6)  both  the  nonlinear 
surface  shapes  and  the  predicted  peak  and  trough  growth  rates  are  are  identical.  In  addition,  convergence 
studies,  carried  out  on  grids  of  various  sizes,  indicate  a  second-order  convergence  rate  consistent  with  the 
linear  elements  used  in  the  BEM  implementation. 

A  comparison  has  also  been  generated  for  highly-distorted  surface  shape  measurements  of  Rutland  and 
Jameson.14  Figure  4  provides  the  comparison  between  the  experimental  wave  shape  reported  by  Rutland 
and  Jameson  and  that  calculated  by  the  BEM.  Two  different  density  ratios  (c  =  0;  0.00129)  are  examined 
for  We  =  109.6;  half  a  wavelength  of  the  surface  is  plotted.  Calculation  of  the  volume  of  fluid  present  in  the 
Rutland  and  Jameson  result  indicate  a  slightly  lower  wave  number  ( k  =  0.237)  than  that  quoted  ( k  =  0.25). 
This  difference  could  easily  be  attributed  to  “accordion  effects”  of  wave  stretching  due  to  the  fact  that  the 
domain  is  in  fact  finite.  Using  k  —  0.237,  results  in  Fig.  4  were  obtained  by  matching  the  height  of  the 
crest  and  plotting  the  resulting  waveform.  The  results  indicate  that  even  at  the  relatively  low  We  value, 
aerodynamic  interactions  with  the  gas  phase  are  important  in  determining  the  nonlinear  wave  shapes.  The 
“swelling”  in  the  trough  region  of  the  experimental  data  is  also  predicted  for  the  e  —  0.00129  case  using 
the  present  model,  although  not  to  the  extent  measured  by  Rutland  and  Jameson.  Nevertheless,  the  result 
including  aerodynamic  effects  is  clearly  superior  to  the  calculation  ignoring  gas  presence.  Errors  between  the 
e  =  0.00129  result  and  the  experiment  can  be  attributed  to  the  nonperiodicity  of  the  actual  measurement  as 
well  as  viscous  effects. 

IV.  Results  and  Discussion 

Results  of  calculations  using  the  model  are  shown  in  Figs.  5-11.  Figure  5  presents  a  typical  growth 
rate  for  the  surface.  In  this  figure  the  square  of  the  growth  rate  (a;)  is  plotted  versus  the  nondimensional 
time  for  two  locations  on  the  jet  surface  for  We  =  850,  e  =  0.00129,  k  =  1.07,  and  89  nodes  along  the  jet 
surface.  The  growth  rate  is  plotted  for  both  the  initial  peak  ( kz  =  0)  and  the  initial  trough  ( kz  =  7r)  of  the 
disturbance  along  with  the  value  predicted  by  the  linear  theory.  The  computed  growth  rates  for  both  the 
peak  and  the  trough  agree  very  well  with  the  linear  theory  for  about  the  first  250  time  increments.  At  this 
point,  the  growth  of  the  jet  departs  from  the  linear  theory  and  becomes  nonlinear.  The  growth  rate  of  the 
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trough  dramatically  decreases  once  the  jet  is  in  the  nonlinear  portion  of  the  process,  and  it  is  this  decrease 
in  the  growth  rate  of  the  trough  (equivalent  to  a  decrease  in  the  radial  velocity  of  the  trough)  that  leads  to 
the  formation  of  the  main  and  the  satellite  drops.  If  the  trough  did  not  undergo  this  decrease  in  the  growth 
rate,  the  jet  would  form  only  one  droplet  per  wavelength  when  it  pinched.  The  “wiggles”  on  the  growth 
rates  during  the  first  50  time  increments  are  a  result  of  the  method  of  treating  the  corner  nodes  in  the  gas 
domain.  The  magnitude  of  these  oscillations  is  reduced  by  increasing  the  number  of  nodes  along  the  wave, 
or  decreasing  either  the  initial  time  step  the  initial  disturbance  amplitude.  The  presence  of  these  oscillations 
has  a  negligible  impact  on  either  the  breakup  length  of  the  jet  or  on  the  drop  sizes  formed. 

Figure  6  presents  the  nonlinear  jet  evolution  in  the  first  wind-induced  regime.  This  plot  is  for  the  same 
conditions  as  the  growth  rate  plot  (Fig.  5).  In  this  figure,  the  surface  shape  is  given  at  three  different  time 
steps  during  the  evolution,  with  two  complete  waves  of  the  disturbance  shown.  The  first  surface  shape  is 
at  t  =  250,  and  corresponds  to  a  point  just  after  the  jet  enters  into  the  nonlinear  portion  of  its  growth.  At 
this  point,  the  maximum  and  the  minimum  radii  of  the  surface  are  at  the  points  that  correspond  to  the 
initial  peaks  and  troughs  of  the  wave.  The  second  surface  shape  is  given  at  t  =  300,  well  into  the  nonlinear 
portion  of  the  jet  growth.  At  this  time,  the  minimum  radius  point  on  the  surface  no  longer  corresponds  to 
the  initial  trough.  The  trough  area  has  flattened  out  and  there  are  now  two  points  of  minimum  radius  per 
wavelength,  one  at  each  end  of  the  flattened  trough  area.  The  final  surface  shape  shown  is  for  t  =  312,  just 
before  the  pinching  process  occurs.  At  this  time,  the  larger  main  drops  are  separated  by  thin  ligaments  of 
fluid  that  make  up  the  satellite  drops.  The  swelling  phenomena  is  not  as  pronounced  as  that  in  Fig.  4  for 
this  relatively  short  wavelength  disturbance. 

Using  this  model,  we  have  identified  the  region  in  which  jet  breakup  transitions  from  the  first  to  second 
wind-induced  regimes.  Figure  7  gives  the  nonlinear  jet  evolution  for  the  transition  region  for  We  =  1800, 
€  =  0.00129,  k  =  1.68,  and  79  nodes  on  the  interface.  Again,  the  jet  surface  is  shown  at  three  different 
times  during  its  evolution.  The  evolution  of  the  trough  through  the  times  shown  is  similar  to  the  trough  in 
the  first  wind-induced  regime.  The  evolution  of  the  peak  however,  is  significantly  different  than  it  is  for  the 
first  wind-induced  case.  At  t  —  200,  the  peak  is  still  rounded.  As  the  surface  evolves  to  t  —  225,  the  peak 
loses  its  rounded  shape  and  takes  on  more  of  a  “peaked”  appearance.  Between  t  =  225  and  t  =  241,  the 
peak  undergoes  even  a  more  dramatic  change  in  appearance.  The  peak  loses  all  of  the  rounded  shape  and  it 
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now  forms  a  “ring”  type  structure  around  the  edge  of  the  jet.  It  is  the  formation  of  this  ring  structure  that 
leads  to  droplet  shedding  from  the  periphery  of  the  jet  in  the  second  wind-induced  region.  When  the  jet 
pinches  in  this  transition  regime,  the  satellite  drops  still  have  the  same  appearance  as  they  did  in  the  first 
wind-induced  regime,  while  the  main  drops  are  no  longer  spherical  and  are  drawn  out  into  a  “ring”  around 
the  edge  of  the  jet. 

The  nonlinear  jet  evolution  in  the  second  wind-induced  regime  is  given  in  Figure  8  for  We  =  5000, 
e  =  0.00129,  k  =  4.23,  and  89  nodes  on  the  interface.  In  this  figure,  the  jet  surface  shape  is  given  at  four 
different  times  during  its  evolution.  The  evolution  of  the  surface  shape  is  significantly  different  than  it  was 
for  the  lower  Weber  number  jets.  At  the  first  time  shown  (<  —  50),  the  jet  has  undergone  almost  no  motion 
of  the  surface  and  is  still  sinusoidal.  At  t  =  80,  the  surface  still  appears  sinusoidal  although  the  peak  has 
started  pulling  outward.  At  t  —  90,  the  peak  has  acquired  a  spiked  shape  and  is  moving  outward  with  a 
rapid  velocity.  When  the  jet  pinches  at  t  —  100,  the  peak  has  pulled  up  a  significant  distance  from  the 
surface  and  has  a  definite  “ring”  structure  at  the  edge.  When  this  ring  of  fluid  is  shed  from  the  periphery 
of  the  jet,  the  trough  is  not  close  to  the  centerline  of  the  jet.  Once  the  ring  of  fluid  is  shed,  it  will  further 
break  up  into  a  series  of  smaller  drops  around  the  edge  of  the  jet,  and  these  small  drops  are  characteristic 
of  the  second  wind-induced  regime. 

At  this  point,  it  is  prudent  to  caution  the  reader  of  the  effects  of  gas-phase  viscosity  on  the  predictions 
in  Figs.  6-8.  The  most  important  real  fluid  effect  is  the  flow  separation  which  will  undoubtedly  occur  in  the 
highly-distorted  shapes  obtained  at  large  times.  The  net  drag  force  induced  by  this  flow  separation  (Fq) 
can  be  expressed  in  terms  of  the  dynamic  pressure  and  the  height  of  the  wave: 

Fo  oc  pgU2tj  (22) 

whereas,  the  net  force  generated  due  to  surface  tension  at  the  interface  is: 

Fa  <x  <r  (23) 

so  that  viscous  effects  will  become  important  as  the  ratio: 

Fd/F„  m  eWe—  (24) 

a 

approaches  unity.  This  condition  is  roughly  attained  at  times  t=300,  200,  and  80,  respectively,  in  Figs. 
6-8.  Beyond  these  times,  we  expect  a  noticible  deflection  of  the  peaks  in  the  downstream  direction  for  the 
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case  of  a  viscous  fluid.  However,  we  note  that  gas  separation  is  not  an  essential  mechanism  to  explaining 
the  gas/liquid  interaction,  particularly  at  conditions  where  the  ratio  Fq/F a  is  small  compared  to  unity. 
A  complete  viscous  analysis  would  be  quite  demanding  computationally,  since  we  could  no  longer  employ 
periodic  boundary  conditions  due  to  differences  in  boundary  layer  thickness  at  either  end  of  a  given  wave. 

Figure  9  compares  the  drop  sizes  predicted  by  the  model  with  those  measured  by  Rutland  and  Jameson8 
and  Lafrance7.  Unfortunately,  all  of  these  data  lie  well  within  the  Rayleigh  regime;  we  were  unable  to 
find  analogous  measurements  at  higher  jet  velocities.  Our  calculations  slightly  underpredict  the  size  of  the 
main  drops,  resulting  in  a  corresponding  overprediction  of  the  satellite  drop  sizes.  This  result  is  essentially 
the  same  as  that  calculated  by  Mansour  and  Lundgren  (with  the  exception  that  actual  gas  density  and  jet 
velocity  were  used  for  each  datapoint)  due  to  the  fact  that  low  speed  jets  were  utilized  in  generating  the 
data.  There  were  some  subtle  differences  obtained  by  including  gas  phase  interactions,  but  in  general  these 
differences  were  of  the  order  of  1-2%  in  droplet  diameters. 

Figure  10  depicts  the  radius  of  the  main  and  the  satellite  drops  versus  the  wave  number  of  the  disturbance 
for  e  =  0.00129  and  We  =  2,  850,  and  1500.  For  We  =  2,  the  main  and  satellite  drops  axe  of  equal  size 
for  k  =  0.28.  Below  this  wave  number,  the  satellite  drops  are  larger  than  the  main  drops,  while  above  this 
wavenumber,  the  main  drops  are  larger  than  the  satellite  drops.  At  increased  jet  velocities  (i.e.  Weber 
numbers),  the  main  drop  radius  decreases  and  the  satellite  drop  correspondingly  increases  in  size.  This 
effect  is  attributed  to  the  importance  of  the  swelling  phenomena  described  in  Fig.  4.  At  higher  We  values, 
the  swelling  in  the  trough  region  is  more  dramatic  (particularly  for  lower  k  values),  thus  moving  the  pinch 
location  closer  to  the  peak  and  increasing  the  size  of  the  satellite  droplet. 

A  series  of  calculations  were  performed  at  various  density  ratios  in  order  to  identify  the  the  transition 
point  between  the  flow  regimes.  The  transition  between  the  Rayleigh  and  the  first  wind-induced  regimes 
occurs  at  the  longest  breakup  length  jet.  Using  a  linear  analysis,  one  can  show  that  aerodynamic  forces  are 
of  the  same  order  as  the  capillary  force  when  the  gas-based  Weber  number  (Wetf)  is  equal  to  unity.  This 
condition  corresponds  to  the  transition  between  Rayleigh  and  first  wind-induced  regimes.  The  transition 
between  first  and  second  wind-induced  regimes  is  only  evident  under  large  deformations  and  can  only  be 
predicted  using  a  complete  nonlinear  analysis.  This  transition  is  important  since  it  signals  the  point  at  which 
atomization  occurs  at  the  jet  periphery,  rather  than  at  the  centerline. 
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Figure  11  shows  results  of  the  transitions  studies  in  terms  of  a  the  jet  breakup  length  versus  We^  for 
various  c  values.  For  each  density  ratio,  the  transition  point  between  the  Rayleigh  and  the  first  wind-induced 
regimes  occurs  at  We^  =  1.0,  as  predicted  by  linear  theory.  This  point  corresponds  to  the  maximum  breakup 
length  of  the  jet  as  depicted  in  Fig.  11.  Agreement  with  the  linear  result  We^  =  1  is  not  surprising,  since 
the  growth  rate  is  near  the  linear  result  for  the  bulk  of  the  process  (e.g.  Fig.  5). 

For  a  given  density  ratio,  Wey  was  increased  to  the  point  where  atomization  began  to  occur  at  the 
periphery  of  the  jet.  This  point,  noted  by  the  termination  of  each  of  the  curves  in  Fig.  11,  signals  transition 
from  the  first  to  the  second  wind-induced  regimes.  Results  indicate  that  this  transition  occurs  at  Wetf  «  2.5 
over  the  large  range  of  density  ratios  studied.  Figure  11  also  shows  the  effect  of  increasing  the  density  ratio 
on  the  breakup  length  of  the  jet.  As  density  ratio  increases,  the  aerodynamic  forces  become  more  important, 
and  the  jet  becomes  more  unstable  and  has  a  shorter  length  upon  breakup  as  suggested  by  the  linear  theory. 
VII.  Conclusions 

Nonlinear  atomization  processes  have  been  modeled  for  a  liquid  jet  under  the  influence  of  both  capillary 
forces  and  aerodynamic  interactions  with  an  external  gas.  The  model  is  capable  of  predicting  the  nonlinear 
evolution  of  the  gas/liquid  interface  within  both  wind-induced  regimes.  Validations  against  linear  theory 
suggest  adequate  agreement  in  the  case  where  initial  jet  deflections  are  assumed  to  be  0.1%  of  the  initial 
jet  radius.  Significant  nonlinear  effects  are  present  for  initial  deflections  as  small  as  1%  of  the  jet  radius. 
Comparison  with  highly  distorted  surface  shape  measurements  indicate  that  the  gas  phase  interaction  is 
important  even  at  relatively  low  jet  velocities.  The  presence  of  the  gas  leads  to  a  “swelling”  in  the  trough 
region  of  the  wave  which  is  predicted  by  the  current  model.  Aerodynamic  interactions  had  very  little  effect 
on  predicted  droplet  sizes  for  low  speed  jets  within  the  Rayleigh  regime.  At  higher  velocities,  a  decrease  in 
main  drop  size  (with  an  attendant  increase  in  satellite  drop  size)  is  predicted  by  the  model.  This  behavior 
is  attributed  to  the  swelling  phenomena  which  effectively  drives  the  pinch  location  toward  the  main  droplet 
at  higher  jet  velocities.  Results  also  indicate  that  transitions  from  Rayleigh  to  first  wind-induced  and  from 
first  to  second  wind  induced  regimes  occur  at  We^  =  1  and  We^  «  2.5,  respectively. 
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Figure  Captions 

FIG.  1.  Regimes  of  Jet  Breakup:  a)  Rayleigh,  b)  First  Wind-Induced,  c)  Second  Wind-Induced,  d)  Atom¬ 
ization 

FIG.  2.  Schematic  of  Computational  Domain  Denoting  Boundary  Conditions 

FIG.  3.  Comparison  of  BEM  to  Linear  Theory  for  €  =  0.001  and  (a)  We  =  2,  1500  and  (b)  We  =  15,000 
FIG.  4.  Comparison  to  Wave  Shapes  of  Rutland  and  Jameson14,  k  =  0.237,  We  =  109.6,  e  =  0.0,  0.00129 
FIG.  5.  Growth  Rate  for  the  Initial  Peak  and  Trough  in  the  First  Wind-Induced  Regime,  e  =  0.00129, 
r]  =  0.004,  k  =  1.07,  We  =  850 

FIG.  6.  Nonlinear  Jet  Evolution  in  the  First  Wind-Induced  Regime,  e  =  0.00129,  k  =  1.07,  rj  =  0.004,  We 
=  850 

FIG.  7.  Nonlinear  Jet  Evolution  in  the  Transition  Region,  e  =  0.00129,  t]  =  0.01,  k  =  1.68,  We  =  1800 
FIG.  8.  Nonlinear  Jet  Evolution  in  the  Second  Wind-Induced  Regime,  e  —  0.00129,  r?  =  0.01,  k  —  4.23,  We 
=  5000 

FIG.  9.  Comparison  of  the  BEM  to  the  Drop  Sizes  of  Rutland  and  Jameson8  (o)  and  Lafrance7  (x)  for  Main 
and  Satellite  Drops,  e  =  0.00129 

FIG.  10.  Drop  Size  versus  Wave  Number  for  e  =  0.00129,  We  =  2,  850,  1500 

FIG.  11.  Breakup  Length  versus  Gas  Based  Weber  Number  for  e  =  0.00129  ( — ),  0.00165  ( — ),  0.005  ( — ), 
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Abstract 

A  Boundary  Element  Method  (BEM)  has  been  developed  to  investigate  the  nonlinear  evolution  of  the  surface 
of  liquid  jets  injected  from  circular  orifices.  The  current  treatment  focuses  on  the  low  speed  jet  in  which 
liquid  viscosity  and  gas-phase  pressure  interactions  are  of  minor  importance.  The  model  properly  reflects 
both  the  presence  of  the  orifice  as  well  as  the  internal  flow  geometry.  The  model  has  been  validated  against 
nonlinear  droplet  oscillation  calculations  of  other  researchers.  Results  are  presented  for  flows  within  the 
dripping  and  Rayleigh  breakup  regimes. 

Introduction 

The  atomization  of  liquid  jets  is  a  fundamental  problem  of  fluid  mechanics.  The  behavior  of  this  complex 
process  is  of  engineering  interest  in  numerous  applications  such  as  fuel  injectors,  paint  applicators,  and  ink 
jet  printers  (to  name  just  a  few).  For  this  reason,  the  atomization  process  has  been  studied  by  a  large  number 
of  researchers,  both  experimentally  and  analytically. 

For  those  seeking  analytic  solutions  related  to  this  problem,  the  nonlinearity  associated  with  the  free 
surface  boundary  condition  represents  a  formidible  obstacle.  Many  early  works  effectively  eliminated  this 
nonlinearity  by  implementing  a  linear  stability  analysis  on  an  infinite  cylindrical  column  of  fluid  in  which  free 
surface  boundary  conditions  are  employed  on  the  undisturbed  interface.  A  good  summary  of  these  works 
is  provided  in  Reitz  and  Bracco  [1]  and  Sterling  and  Sleicher  [2].  More  recently,  Leib  and  Goldstein  [3] 
effectively  included  the  presence  of  the  orifice  in  a  linear  analysis  which  also  included  the  effect  of  injection 
velocity  profile  on  jet  stability.  Computational  models  [4,5]  have  also  made  use  of  results  from  the  linear 
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analyses,  combined  with  one-dimensional  mass  and  momentum  balances  to  effectively  model  the  atomization 


process. 

The  increased  computational  resources  available  during  recent  years  has  led  analysts  to  the  development 
of  nonlinear  models  of  these  free  surface  flows.  Traditional  computational  fluid  dynamics  (CFD)  approaches 
are  difficult  to  implement  in  these  problems  due  to  the  large  distortion  of  the  grid  which  occurs  near  the 
droplet  pinch-off  event.  Deng  and  Jeng  [6]  have  had  success  in  modeling  a  liquid  droplet  (in  a  gas  flow)  in 
this  manner,  but  the  solution  accuracy  is  greatly  diminished  if  the  droplet  surface  is  highly  distorted. 

Other  researchers  have  employed  the  Volume  of  Fluid  (VOF)  technique  pioneered  by  Hirt  and  Nichols  [7]. 
This  technique  has  the  advantage  of  being  able  to  treat  a  large  number  of  droplets  which  are  interpolated 
from  a  fixed  (generally  cartesian)  grid.  However,  the  free  surface  must  be  interpolated  from  existing  nodal 
locations  which  leads  to  great  difficulties  in  obtaining  accurate  representation  of  surface  curvature.  While 
this  approach  is  useful  for  modeling  a  droplet  field  within  a  spray  [8,9],  early  stages  of  the  atomization 
process  are  more  difficult  to  handle. 

Recently,  Boundary  Element  Methods  (BEMs)  have  been  utilized  to  accurately  model  free  surface  flows 
involving  droplets  [10]  and  liquid  jets  [11-13]  for  free  surface  flows.  Since  these  techniques  utilize  a  “grid” 
of  nodes  placed  only  on  the  boundary  of  the  domain,  they  can  handle  large  surface  distortions  and  droplet 
pinchoff  without  complications.  Previous  liquid  jet  studies  have  revealed  the  presence  of  satellite  droplets  in 
the  Rayleigh  [11]  and  first  wind-induced  [12]  regimes  through  the  examination  of  a  single  wave  in  a  liquid 
jet  of  infinite  axial  extent.  However,  these  models  have  not  investigated  the  effects  of  the  presence  of  the 
orifice  due  to  the  infinite  jet  assumption.  In  very  low  speed  (dripping)  and  high  speed  (atomization)  flows, 
the  orifice  interaction  can  be  important  since  breakup  occurs  in  close  proximity  to  the  exit  plane. 

In  this  paper,  we  describe  a  BEM  capable  of  addressing  the  finite- length  liquid  jet.  The  model  is  described 
in  the  following  section.  Following  this  description,  we  discuss  validation  of  the  code  for  a  nonlinear  droplet 
oscillation  problem.  Nonlinear  results  are  discussed  for  dripping  flows  as  well  as  flow  within  the  Rayleigh 
regime.  Finally,  a  model  of  a  fountain  is  briefly  described. 
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Model  Development 


We  assume  an  axisymmetric,  incompressible,  inviscid  flow  in  which  gas  pressure  variations  are  negligible. 
The  flow  regime  under  which  these  assumptions  are  valid  is  one  in  which  the  liquid  velocity  is  low  since 
viscous  forces  will  become  important  when  the  ratio  fiv/a  approaches  unity.  Under  these  assumptions,  the 
unsteady  liquid  flow  is  described  by  Laplace’s  equation: 

V2«A  =  0  (1) 


where  <j>  is  the  velocity  potential  defined  as  the  function  whose  gradient  is  simply  the  velocity,  i.e.  V  •  <i>  =  v . 
Following  Liggett  and  Liu  [14],  the  BEM  formulation  of  Eq.  (1)  becomes: 


+  J[<1>^  -  qG}dr  =  o  (2) 

where  <£(rj)  is  the  value  of  the  potential  at  a  point  rj,  T  denotes  the  boundary  of  the  domain,  and  G  is  the 
Green’s  function  corresponding  to  Laplace’s  equation.  Since  Eq.  (2)  involves  an  integration  only  around  the 
boundary,  we  need  not  discretize  the  entire  domain.  It  is  presumed  that  either  <f>  or  q  =  d<j>/dn  is  specified 
at  each  “node”  on  the  boundary.  Here  n  is  the  outward  normal  to  the  boundary  so  that  q  represents  the 
velocity  normal  to  the  boundary.  The  quantity  a  in  Eq.  (2)  results  from  singularities  introduced  as  the 
integration  passes  over  the  boundary  point,  r\.  For  an  axisymmetric  formulation,  a  is  equal  to  the  angle 
between  successive  nodes  on  the  boundary. 

The  geometry  for  the  integration  process  is  shown  in  Fig.  ??.  If  we  let  r  and  z  denote  radial  and  axial 
coordinates,  respectively,  and  denote  the  base  point  with  subscript  “i”,  the  Green’s  function  solution  to  the 
axisymmetric  Laplacian  can  be  written: 


q  _  4r/iT(p) 

V5 

where 

a  =  (r  +  r,)2  +  (z  -  z{)2 


(3) 

(4) 


and 


a 

and  K(p)  is  the  complete  elliptic  integral  of  the  first  kind  (for  efficiency  this  quantity  is  calculated  using  a 
curve  fit  [15]. 
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Physically,  the  Green’s  function  represents  the  contribution  of  a  source  (at  location  r)  of  strength  2x  as 
sensed  at  the  “base  point”,  rj.  While  the  quantity  dG/dn  represents  the  contribution  of  a  doublet  as  sensed 
at  the  base  point.  By  properly  weighting  all  these  source  and  doublet  strengths,  we  can  solve  Laplace’s 
equation  for  arbitrary  boundary  conditions.  Derivation  of  dG/dn  for  this  axisymmetric  problem  is  quite 
involved.  Gipson  [16]  gives: 


dG 

dn 


-1 

2  y/a. 


nrK(p)  + 


Ejp) 

(r  ~  r<)2  +  (z  —  2j)2 


[((r  ~  r»)2  ~  (z  ~  zi )2)  +  2rhl(z  - 


(6) 


with  the  quantities  rir  and  hz  being  the  components  of  the  outward  normal  unit  vector  in  the  radial  and 
axial  directions,  respectively,  while  E(p)  is  the  complete  elliptic  integral  of  the  second  kind. 

Performing  the  desired  integration  around  the  boundary  of  the  domain  results  in  an  algebraic  system  of 
equations.  In  order  to  perform  this  integration,  it  is  necessary  to  assume  a  behavior  for  <£  and  q  over  the 
length  of  an  element.  We  assume  that  both  <j>  and  q  vary  linearly  over  the  length  of  an  element.  With  this 
assumption,  the  values  for  the  potential  of  the  element  bounded  by  nodes  j  and  j  +  1  is: 


4>  =  <t>j  ■ 


r,+i  -  r 


+  4>j+i 


r-r, 


(7) 


Tj+i-Tj  7  ri+1-ri 

with  an  analogous  expression  for  the  velocity,  q.  Using  linear  elements,  Eq.  (2)  is  divided  into  two  separate 
integrals: 

r  n  —  1 

/  qGdT  =  £  [qjh(iJ)  +  qj+ih(i,j)}  (8) 

Jr 


and 


jr  =  £  fahiij)  +  4>J+iU(i,j)} 


(9) 


where 


and 


h 


h 


h  = 


h  = 


rJ+i  -r 

ri+i  -  r; 


GdT 


GdT 


r-r, 
r;+l  -  r, 


rj+1  -  r  dG 
r,+i  —  r,  dn 


dr 


r  -  r,  dG 
r,+i  -  r;  dn 


dr 


(10) 

(11) 

(12) 


(13) 
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The  integrals  I\  through  /4  are  calculated  using  four  point  Gaussian  Quadrature  [15].  Using  the  integrals 
through  /4,  we  create  D  and  S  matrices  which  contain  doublet  and  source  contributions,  respectively.  From 
Eqs.  (10)  and  (11),  we  can  write: 

SiJ  =  11  i.i  +  DiJ  =  h'J  +  (14) 

The  integration  of  Eq.  (2)  results  in  a  system  of  algebraic  equations  of  the  form: 

{Dm  =  [S]{q}  (15) 

with  the  values  of  a  incorporated  into  diagonal  elements  of  the  D  matrix.  To  solve  this  system  of  equations, 
all  the  unknown  terms  are  “pivoted”  by  interchanging  the  columns  of  the  appropriate  matrices.  With  all 
the  known  values  on  the  right  hand  side,  the  algebraic  system  of  equations  can  be  written  in  standard  form: 

[A]{x]  =  {6}  (16) 

where  x  is  the  vector  of  unknown  boundary  values  of  <j>  or  q .  The  resulting  matrix  A  is  fully  populated. 
We  solve  this  system  of  equations  by  a  Crout  Method  LU  decomposition  taken  from  Numerical  Recipes  in 
FORTRAN  [17]. 

This  method  requires  n3 / 3  operations  to  solve  the  system,  while  the  popular  Gaussian  elimination  method 
requires  2n3/3  operations.  This  reduction  in  the  number  of  operations  by  a  factor  of  two  is  a  significant 
improvement  in  the  computation  time  required  for  large  systems  of  equations  such  as  those  solved  here. 

Centerline,  Corner,  and  Self-Adjoint  Nodes 

The  calculation  of  the  integrals  in  Eqs.  (10)-(13)  is  straightforward  except  for  three  special  cases:  nodes  on 
the  centerline  of  the  domain,  nodes  on  the  corners  of  the  domain,  and  self-adjoint  nodes. 

^From  the  Green’s  function,  Eq.  (3),  both  G  and  dG/dn  vanish  along  the  centerline  (r  =  0)  of  the 
domain.  Therefore,  no  nodes  need  to  be  placed  along  the  centerline  of  the  domain  except  at  the  ends.  At 
the  endpoints,  the  values  of  G  and  dG/dn  are  not  zero  since  the  integration  at  these  points  is  not  directly 
along  the  centerline.  At  the  first  node,  the  7i(l,  1)  and  ^(1,1)  terms  are  equal  to  7rr2,  where  r2  is  the  radial 
distance  to  the  second  node.  The  ^(l,  1)  and  J4(l,  1)  terms  are  zero.  An  examination  of  the  last  node  in 
the  domain  yields  an  equivalent  result. 
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At  the  corners  of  the  domain,  the  unit  normal  and  normal  velocity  are  multivalued.  In  Fig.  ??,  for 
example,  the  node  labeled  k  has  both  a  horizontal  and  vertical  velocity.  The  boundary  element  method 
only  allows  for  the  calculation  of  one  velocity  or  potential  at  any  one  node.  This  characteristic  of  the  BEM 
requires  special  treatment  of  corner  nodes.  At  this  corner  node,  ifc,  the  general  S  matrix  calculation  can  be 
written: 


Sq  =...  +  [h(i,k-  1)  +  -  2)]gt_! 

+  [I2(i,k-l)]q;  +  [h{i,k)]qt  (17) 

+  [-M*>^+l)  +  fc)]9t+i  + 


where  k  —  1  and  k  +  1  are  the  node  before  and  the  node  after  the  corner,  respectively.  The  elements  of  the 
S  matrix  are  modified  at  the  corners  by  using  values  of  I\  and  1 2  consistent  with  the  desired  solution.  For 
corners  where  the  velocity  is  zero,  such  as  an  orifice  exit  or  nozzle  inflow,  the  element  of  the  S  matrix  is 
equal  to  the  appropriate  value  of  I2  only. 

Self-adjoint  nodes  are  defined  as  the  nodes  in  the  integration  where  i  =  j  and  i  =  j  - f  1.  In  these  two 
cases,  the  elliptic  integral  K(p)  is  singular  due  to  p  having  a  value  of  zero  at  one  end  of  the  integration.  To 
treat  this  singularity,  the  curve  fit  used  for  K(p)  from  Abramowitz  and  Stegun  [15]  is  divided  into  singular 
and  nonsingular  parts: 


K(p)  -  K{p)*  -  b0  ln(p)  =  K(pY  -  1  In 


•(r-r,-)2  +  (z-;Q2- 

>  +  r,-)2  +  (*-z,-)2. 


(18) 


where  K{p)*  represents  the  nonsingular  part  of  the  elliptic  integral.  Substituting  into  Eq.  (10),  the  new 
expression  can  be  divided  into  singular  and  nonsingular  parts.  The  nonsingular  parts  are  treated  similarly  to 
the  nonsingular  equation  for  I\.  For  the  singular  parts,  we  make  the  substitution  f  =  (r  —  rj)/AFj,  which 
results  in  a  third  equation  that  can  also  be  divided  into  singular  and  nonsingular  parts.  The  singular  part  of 
this  new  equation  is  computed  with  special  Gaussian  quadrature  [15].  The  final  result  at  the  point  i  =  j  is: 

Ii(i,  i)  =  -  Ari  J\l  -  ln^  (19) 


where: 

(G),  =  il(w+1n^) 


(20) 
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The  same  approach  is  applied  to  the  I2  through  /4  terms,  and  yields: 


h(i,  0  =  (G).  ^  ~  Ar;  fo  (21) 

h{i' f)  =  ^  nr  (Sr) ,  ■* + ar-  i 1(1  '"<«  <22> 

u(i'  ° = ^  L  V  ® « + Ar*  i‘  ^  ‘23> 

1*  +  ln-&ll 

(24) 

The  self-adjoint  terms  at  i  =  j  -f  1  are  examined  with  the  same  procedure,  and  results  similar  to  Eqs.  (19) 


where, 


^  — _ -E(p) 

>-r.-)2-(z- 

z<)2]  hr  +  2n,r(z  -  z.)  +  hr 

(• K(P )*  +  ln  aT") 

),  -  v*E{p) 

(r-r.)2+(*-z,)2 

and  (21)-(23)  are  obtained. 


Free  Surface  Module 

For  nodes  lying  outside  the  orifice,  boundary  conditions  for  the  “Laplace  Solver”  discussed  above  are  derived 
from  an  examination  of  the  free  surface  boundary  condition.  This  condition  will  enable  us  to  determine  the 
time-dependent  variation  of  <f>  on  the  free  surface.  We  employ  a  procedure  similar  to  that  of  Longuet-Higgins 
and  Cokelet  [18]  to  effectively  track  the  free  surface.  If  we  permit  surface  nodes  to  move  at  the  instantaneous 
velocity,  flow  kinematics  require: 

D1_d±  Dr^dt 

Dt~  dz  Dt  dr  1 } 

where  the  notation  DQ/Dt  implies  a  material  or  Lagrangian  derivative.  Recognizing  that  our  BEM  solver 
will  return  velocities  normal  to  the  surface,  we  employ  the  velocity  transformations: 

It  =  ltsin^  +  qcos(^  f“  =  “  ?s*n(/?)  (26) 

where  /?  is  the  local  wave  slope  and  d<f>/ds  is  the  velocity  tangential  to  the  local  surface,  as  shown  in  Fig. 
??.  This  tangential  velocity  is  calculated  using  5-point  centered  differences  on  <£,  except  for  nodes  adjacent 
to  the  orifice,  where  a  3-point  formula  is  employed. 

The  local  wave  slope,  /?,  is  calculated  following  the  formulation  of  Medina  [19].  For  each  node,  a  parabola 
is  defined  such  that  it  passes  through  the  previous  node,  the  node  in  question,  and  the  following  node.  The 
slope  of  the  surface  is  given  by  the  tangent  to  the  parabola  at  the  central  node. 
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We  obtain  a  means  of  updating  <j>  on  the  surface  boundary  by  considering  the  unsteady  Bernoulli  relation. 
For  our  problem,  after  transforming  the  time  derivative  of  <j>  from  Eulerian  to  Lagrangian  representation, 
the  dimensionless  representation  of  this  equation  is: 


D<f>  l/t_.x2  K  B° 

Dt  2K  We  We 


(27) 


where  We  =  ( pv2a)/a  is  the  Weber  number,  Bo  —  (pga2)/a  is  the  Bond  number,  and  k  represents  the  local 
surface  curvature.  Here  p  is  fluid  density,  a  is  the  radius  of  the  orifice,  a  is  the  fluid  surface  tension,  and  g 
is  the  acceleration  due  to  gravity. 

The  curvature  of  the  surface  is  calculated  using  a  parametric  form  from  Smirnov  [20]. 


k  = 


\/(r,)2  +  (*')2  [(r'f  +  (z')2] 1 


(28) 


where  the  derivative  terms  are  with  respect  to  the  distance  along  the  surface.  Note  that  (r')2  -f  (z')2  =  1 
for  derivatives  with  respect  to  the  distance  along  the  surface.  We  calculate  these  derivatives  using  five  point 
centered  difference  approximations  for  nonuniform  grid  spacing. 


Integration 

In  order  to  integrate  Eqs.  (25)  and  (27),  we  employ  a  fourth-order  Runge-Kutta  method  [21].  The  accuracy  of 
this  method  requires  four  evaluations  of  the  function  for  each  step  in  time.  Even  with  this  large  computational 
cost,  the  method  has  a  significant  benefit  to  the  problem  solved  here.  This  benefit  is  in  “self-starting”.  That 
is,  the  method  does  not  require  information  from  previous  time  steps.  This  self-starting  ability  allows  for  the 
addition  or  subtraction  of  nodes,  the  use  of  a  variable  time  step,  and  for  the  regridding  of  the  surface  over 
time.  Other  popular  fourth  order  methods,  such  as  Adams-Bashforth-Moulton,  are  not  self-starting  and  do 
not  allow  for  these  variations  to  occur. 

Medina  [19]  recognized  that  the  selection  of  a  constant  time  step  depends  on  the  most  demanding  con¬ 
ditions  during  the  surface  evolution.  He  found  the  time  when  the  droplets  are  about  to  be  formed  the  most 
restrictive.  We  employ  a  method  to  dynamically  modify  the  time  step  based  on  the  surface  velocities,  so 
that  the  time  step  is  decreased  if  the  velocities  become  too  large.  In  our  method,  no  node  is  allowed  to  move 
more  than  a  specified  fraction  of  the  grid  spacing  every  time  step.  We  calculate  this  local  time  step  at  every 
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node  by  dividing  the  maximum  allowable  distance  by  the  velocity  (the  maximum  of  the  axial  and  radial 
velocities)  at  that  node.  The  time  step  used  is  the  minimum  of  all  calculated  time  steps  and  the  initial  time 
step  set  as  an  input  to  the  program. 

Numerical  Smoothing 

Several  authors  [10,18]  have  noted  the  appearance  of  “zig-zag”  instabilities  for  BEM  calculations  running 
long  periods  of  time.  To  avert  this  nonphysical  kinking  of  the  surface,  we  apply  a  smoothing  function  every 
time  step  which  is  an  adaptation  of  the  function  used  by  Longuet- Higgins  and  Cokelet  [18]: 

—  (1  —  Cs)<j)j  +  —(“^>-2  4-  4- 10^  4-  4 —  <f>j+ 2)  (29) 

where  Cg  is  a  smoothing  constant  between  0  and  1. 

The  value  of  C3  must  be  large  enough  the  smooth  out  the  point-to-point  kinks  in  the  surface,  but  should 
not  add  a  significant  amount  of  damping  to  the  inviscid  solution.  Using  the  simulation  of  nonlinear  droplet 
oscillations,  various  values  of  Cs  were  used  to  determine  the  appropriate  magnitude  of  this  constant.  Values 
as  large  as  0.10  were  found  to  decrease  the  total  energy  of  the  droplet.  Therefore,  we  typically  use  a  value 
near  0.01,  which  provided  good  results  without  damping  the  droplet  energy. 

Regridding 

Since  the  surface  nodes  are  allowed  to  move  with  the  local  surface  velocity,  the  spacing  between  nodes  changes 
in  time.  In  order  to  keep  the  grid  spacing  roughly  equal,  the  surface  is  remeshed  every  time  step  using  cubic 
splines.  Using  the  distance  along  the  surface,  s,  as  the  spline  parameter,  the  spline  coefficients  for  r(s),  *(s), 
and  <f>(s)  are  calculated  by  solving  a  tridiagonal  system  of  equations  for  the  first  derivative  coefficients  [21]. 
The  slope  at  the  ends  of  the  spline  are  calculated  where  appropriate  using  5-point  centered  differences. 

For  the  oscillating  droplet  case,  the  nodes  are  redistributed  keeping  the  total  number  of  nodes  a  constant. 
For  the  liquid  jet,  the  grid  spacing  along  the  jet  is  held  roughly  constant,  allowing  the  number  of  nodes  to 
increase  as  the  jet  issues  from  the  orifice,  and  decrease  as  droplets  are  shed  from  the  calculation.  As  the  liquid 
jet  begins  to  neck  down,  a  node  is  maintained  at  the  point  of  minimum  radius.  For  the  jet  results  presented 
here,  the  grid  spacing  is  kept  near  20%  of  the  orifice  radius.  A  droplet  is  separated  from  the  calculation 


46 


when  the  jet  pinches  below  5%  of  the  orifice  radius.  Numerical  experiments  indicate  that  solutions  are  not 
sensitive  to  either  of  these  tolerance  values. 


Model  Validation  -  Nonlinear  Droplet  Oscillations 

To  validate  the  axisymmetric  BEM  solver  and  free  surface  modules,  we  use  results  for  the  oscillation  of  a 
liquid  droplet  in  the  absence  of  gravity.  Lamb  presents  the  linearized  solution  for  small  oscillations  of  a  drop 
of  liquid  about  the  spherical  form  [22].  His  results  can  be  expressed  (for  a  given  mode,  n): 

r  =  a  +  enP„(cos0)sin(ui°t  +  tj)  (30) 


for  the  surface  shape, 

<t>  -  ~~~  (“)  tnPn(cO80)cOs(w°t  +  T) )  (31) 

for  the  velocity  potential  inside  the  drop,  and 


.0x2  _  n(n  -  l)(n  +  2)<r 
1 n>  ~  pa3 


(32) 


for  the  frequencies.  Here,  cr  is  the  surface  tension,  a  is  the  unperturbed  radius  of  the  drop,  and  Pn(co$0) 
is  the  Legendre  polynomial  of  order  n  with  9  the  polar  angle.  Since  the  oscillating  droplet  does  not  have 
a  characteristic  velocity,  we  normalize  velocities  by  the  quantity  y/<r/{pa).  For  the  second  mode,  the  non- 
dimensional  frequency  is  given  by  =  y/S/We .  Using  25  nodes  on  half  a  droplet  with  no  smoothing  or 
remeshing,  a  radius  of  1  and  a  displacement  of  0.05,  the  measured  frequencies  for  the  second  mode  from  the 
BEM  code  for  Weber  numbers  of  0.5  to  2.0  are  within  one  tenth  of  one  percent  of  the  theoretical  frequency. 

Tsamopoulos  and  Brown  [23]  studied  moderate  oscillations  of  inviscid  droplets  by  expanding  to  second 
order  in  amplitude.  For  the  surface  shape  of  the  n  =  2  mode,  they  found  the  result  for  t  =  0: 


„  1  o/  2  44  144  . 

r_a  +  e2P2+-c2(--  +  — P2  +  — P4)  (33) 

Figure  ??  shows  the  initial  shape  for  an  amplitude  of  62  =  0.40  and  the  shape  returned  by  the  BEM 
calculation  after  half  a  period  of  oscillation. 

Using  45  nodes  along  the  half  droplet,  we  calculated  droplet  oscillation  frequencies  for  the  n  —  2  mode 
for  initial  disturbances  of  €2  =  0.05  to  0.70.  The  results  are  presented  in  Fig.  ??  as  a  percent  frequency 
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shift  from  the  linear  result  versus  aspect  ratio  (length/width)  at  maximum  prolate  shape.  Tsamopoulos  and 
Brown  considered  moderate  amplitude  oscillations  and  only  presented  their  own  data  to  aspect  ratios  under 
2.  In  this  range,  the  BEM  code  returns  frequency  shifts  within  0.5%  of  the  analytic  result.  Extending  the 
analytic  result  to  aspect  ratios  up  to  about  3,  we  find  frequency  shifts  within  about  2%  of  the  expected 
values.  For  these  large  disturbance  oscillations,  however,  higher  order  corrections  are  needed  if  the  droplet 
is  to  oscillate  more  than  a  few  cycles. 

Lundgren  and  Mansour  [10]  studied  the  oscillations  of  drops  with  weak  viscous  effects  using  the  boundary 
integral  method.  They  present  results  for  various  modes  and  initial  conditions  for  inviscid  droplet  oscillations. 
We  find  excellent  comparisons  with  our  results  for  the  same  cases. 


Dripping  Flows 

Figure  ??  shows  a  typical  computational  domain  for  the  finite-length  liquid  jet.  In  this  figure,  we  consider  a 
sharp-edge  orifice  of  some  given  length  and  inside  diameter.  The  model  is  capable  of  handling  axisymmetric 
orifice  geometries  of  arbitrary  dimensions.  We  employ  a  “hybrid”  set  of  boundary  equations:  nodes  inside 
the  orifice  are  fixed  and  have  a  specified  normal  velocity;  nodes  on  the  free  surface  move  with  the  local 
velocity  and  have  a  known  value  of  the  velocity  potential. 

The  initial  shape  of  the  jet  is  assumed  to  be  cylindrical  with  a  spherical  endcap,  moving  with  a  uniform 
axial  velocity.  The  initial  length  of  the  jet  and  the  radius  of  the  endcap  are  inputs  to  the  initialization  routine. 
By  running  the  simulation  through  several  pinches,  subsequent  results  have  been  shown  to  be  insensitive 
to  this  initial  assumed  geometry.  Because  of  the  assumed  initial  geometry,  data  for  the  first  few  droplet 
pinches  is  not  included  in  the  results.  At  low  velocities,  the  breakup  of  a  liquid  jet  occurs  in  the  “dripping 
flow”  regime.  In  this  regime,  the  gas  around  the  jet  can  be  neglected,  surface  tension  and  gravity  are  the 
important  parameters,  and  the  droplet  diameters  are  larger  than  the  orifice  diameter.  Wu  and  Schelly  [24] 
present  experimental  results  on  the  effects  of  surface  tension  and  temperature  on  the  dynamics  of  a  dripping 
faucet.  Their  nozzle  geometry  is  shown  in  Fig.  ??;  the  orifice  is  beveled  45  degrees,  the  inside  length  of  the 
nozzle  is  40  mm,  and  the  inner  and  orifice  radii  are  5  mm  and  1  mm,  respectively.  Using  a  continuously 
variable  flow  rate,  they  recorded  a  “dripping  spectrum”  (drop  interval  as  a  function  of  flow  rate)  for  various 
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solutions  of  water  and  surfactant.  Results  indicate  chaotic  and  bifurcated  (bimodal)  shedding  frequencies 
over  the  range  0.625  <  We  <  1.05  for  the  case  where  Bo  =  0.204.  For  We  >  1.05,  the  droplet  interval 
increases  with  decreasing  flow  rate  which  implies  an  increasing  jet  length  as  We  increases. 

Figure  ??  presents  a  typical  time  history  for  the  droplet  shedding  process.  As  a  droplet  is  shed,  a 
disturbance  is  sent  upstream.  The  interaction  of  this  disturbance  with  the  orifice  has  been  found  to  be 
responsible  for  the  chaotic/bifurcated  behavior  observed  at  low  We  values.  It  is  interesting  to  note  that  a 
deterministic  model  can  explain  chaotic  behavior  of  this  type.  Previous  efforts  [24]  have  focused  on  chaos 
theory  to  explain  the  observed  behavior. 

Figure  ??  shows  the  droplet  sizes  as  a  function  of  the  Weber  number  for  conditions  consistent  with 
Wu  and  Schelly’s  experiment.  We  note  the  chaotic/bifurcation  behavior  over  a  range  0.6  <  We  <  1.4  and 
drop  sizes  between  1.5  and  2.0  jet  radii.  Plotting  a  “dripping  spectrum”  with  drop  interval  as  a  function  of 
Weber  number  as  shown  in  Fig.  ??,  we  find  that  the  chaotic  region  occurs  at  slightly  higher  Weber  numbers 
and  with  shorter  drop  intervals  than  the  results  of  Wu  and  Schelly.  This  difference  is  due  to  the  stabilizing 
influence  of  viscosity  in  the  actual  measurements.  However,  the  general  character  of  the  map  is  reproduced 
using  the  BEM  calculation.  The  minimum  drop  interval  tends  to  be  20  msec  for  all  Weber  numbers  above 
about  We  =  0.6.  The  jet  breakup  length  in  this  range  of  Weber  numbers  is  shown  in  Fig.  ??.  As  We 
(jet  velocity)  is  increased,  the  jet  breakup  point  reaches  a  location  such  that  this  disturbance  has  vanishing 
influence  on  the  orifice  and  the  chaotic/bifurcated  behavior  vanishes.  Above  Weber  numbers  of  1.9,  the 
model  predicts  a  dramatic  increase  in  jet  breakup  length,  signaling  the  start  of  the  Rayleigh  breakup  regime. 


Rayleigh  Breakup 

At  moderate  inflow  velocities,  the  presence  of  the  surrounding  gas  can  still  be  neglected,  and  the  jet  breakup 
length  increases  as  the  inflow  velocity  is  increased.  Experimental  measurements  of  flows  consistent  with 
the  Rayleigh  regime  [25-28]  show  the  presence  of  satellite  drops  formed  during  the  nonlinear  portion  of 
the  atomization  process.  These  experiments  were  performed  by  using  a  forced  oscillation  frequency  thus 
permitting  investigation  of  the  growth  of  instabilities  over  a  range  of  wavelengths. 

To  simulate  these  experimental  results,  we  employ  a  computational  domain  with  fixed  nodes  only  at  the 
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orifice  exit;  the  nozzle  geometry  is  not  modeled.  We  specify  an  inflow  velocity  along  the  entire  orifice  exit  of 
the  form: 

q=zq[l  +  r]sin(kt)]  (34) 

where  q  is  the  mean  inflow  velocity,  17  is  the  amplitude  of  the  oscillatory  disturbance,  and  k  is  the  wave 
number  of  the  perturbation.  The  calculation  begins  with  a  long  column  of  fluid  outside  the  orifice  in  order 
to  insure  that  perturbations  from  the  assumed  spherical  endcap  do  not  influence  the  development  of  waves 
due  to  the  imposed  perturbation. 

Simulations  were  run  for  We  =  100  with  Bo  =  0  and  0.5292.  The  nonzero  Bond  number  is  consistent 
with  water  flow  through  a  2  mm  radius  orifice.  The  inflow  velocity  had  a  substantial  perturbation  of  3%  of 
the  mean  velocity  with  a  wave  number  of  0.7.  We  are  forced  to  consider  perturbations  of  this  magnatude 
due  to  the  large  computational  resources  required  for  smaller  amplitude  simulations;  i.e.  small  perturbations 
lead  to  very  long  jets.  Since  large  jet  lengths  result  in  a  large  number  of  nodes  along  the  surface  and  long 
computation  times,  we  begin  these  simulations  with  a  grid  spacing  of  50%  of  the  orifice  radius.  After  the  jet 
has  shed  the  initial  fluid  outside  the  orifice,  the  grid  spacing  was  reduced  to  20%  of  the  orifice  radius,  and 
the  simulations  were  continued  through  a  number  of  droplet  shedding  events. 

Figure  ??  shows  the  jet  profile  for  Bo  =  0  at  four  equally  spaced  time  intervals  during  the  simulation.  The 
first  profile  occurred  at  t  —  87.0,  and  the  region  near  the  orifice  is  not  shown  in  the  figure.  The  numerical 
simulation  shows  the  formation  of  the  satellite  droplets  found  in  the  experimental  results.  We  calculate 
average  main  and  satellite  droplet  radii  of  1.793  and  0.986,  respectively,  and  a  breakup  length  of  38.8  orifice 
radii.  For  a  wave  number  of  k  =  0.683  and  We  =  112,  Rutland  and  Jameson  [26]  report  droplet  radii  of 
1.75  and  0.58.  The  difference  between  the  BEM  calculation  and  experimental  result  can  be  attributed  to 
the  amplitude  of  the  perturbation.  For  example,  Dressier  [29]  has  shown  that  uniform-sized  drops  can  be 
generated  if  the  perturbation  is  of  sufficient  amplitude  which  implies  that  the  satellite  droplet  grows  in  size 
with  increasing  perturbation  magnatude. 

When  the  effect  of  gravity  is  included  in  the  simulation,  the  disturbance  wavelengths  were  found  to 
increase  along  the  jet.  A  comparison  of  jet  profiles  for  Bo  =  0  and  Bo  =  0.5292  at  t  =  71  is  shown  in 
Figure  ??.  From  the  figure,  it  is  obvious  that  some  stretching  is  occurring  for  substantial  perturbation 
modeled  in  this  case.  We  have  found  that  this  phenomenon  becomes  increasingly  important  as  the  jet 


50 


perturbation  velocity  component  is  increased  and  is  most  pronounced  near  the  orifice.  We  believe  the 
stretching  phenomena  can  be  attributed  to  differences  in  surface  tension  in  the  axial  and  azimuthal  directions. 
The  azimuthal  load  carrying  capability  (roughly  a /a  near  the  orifice)  provides  an  effective  buckling  strength 
for  the  column  of  fluid  as  it  is  accelerated  during  the  positive  portion  of  the  velocity  oscillation.  However, 
the  axial  contribution,  which  would  provide  tensile  strength  for  the  column  during  the  negative  portion  of 
the  velocity  contribution  is  very  small  near  the  orifice  since  curvature  in  this  plane  is  very  nearly  zero.  For 
very  small  perturbations,  this  phenomena  would  not  be  evident  since  the  jet  would  be  under  negliglible  axial 
acceleration  levels. 

We  have  noted  that  this  phenomena  is  present  in  the  measurements  of  Rutland  and  Jameson  [30]  in 
which  they  measured  the  wave  shape  while  driving  a  jet  with  a  fixed  frequency  (and  hence  fixed  wavelength) 
disturbance.  A  5%  increase  in  jet  volume  was  noted  for  a  case  showing  the  wave  profile  just  prior  to  breakup. 
While  some  of  this  discrepancy  may  be  attributed  to  experimental  error,  the  wave  stretching  phenomena 
described  above  could  also  explain  this  measured  result. 

Fountain  Simulation 

By  specifying  a  negative  Bond  number,  gravity  can  be  directed  opposite  to  the  jet  velocity  vector,  producing 
the  simulation  of  a  fountain.  Figure  ??  shows  the  jet  profile  at  equal  time  intervals  for  a  typical  fountain 
simulation.  In  this  case,  we  choose  water  flowing  at  1  m/sec  through  a  1  cm  radius  orifice.  The  nozzle  has 
length/diameter  ratio  of  0.5,  and  an  inlet  corner  rounded  at  20%  of  the  orifice  radius.  Inflow  is  specified 
normal  to  a  3  cm  radius  hemisphere  at  a  rate  to  give  the  desired  outflow  velocity.  As  the  figure  shows,  the 
jet  forms  a  mushroom-shaped  cap  of  fluid  which  is  terminated  in  a  torus. 


Conclusions 

A  numerical  model  using  the  boundary  element  method  has  been  developed  for  the  purpose  of  investigating 
the  nonlinear  evolution  of  low-speed  the  orifice  presence  for  arbitrary  internal  geometry,  and  resolves  the 
free  surface  of  the  jet  with  high  resolves  the  free  surface  of  the  jet  with  high  accuracy  throughout  the 
calculation.  Multiple  droplets  can  be  removed  from  the  simulation  as  they  are  pinched  off  from  the  liquid 
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jet,  and  both  steady  and  unsteady  injection  processes  can  be  handled.  The  mod.  For  dripping  flows,  chaotic 
and  bifurcated  droplet  shedding  is  predicted  using  this  deterministic  model.  For  dripping  flows,  chaotic  and 
bifurcated  droplet  shedding  is  predicted.  We  find  that  this  behavior  is  explained  by  interactions  of  surface 
waves  (generated  during  pinchoff)  with  the  orifice.  Results  indicate  that  an  abrupt  transition  from  dripping 
flow  to  Rayleigh  flow  occurs  at  We  =  1.9.  The  model  predicts  formation  of  iain  and  satellite  droplets  in 
the  Rayleigh  regime  for  oscillating  inflow  velocities.  The  acceleration  due  to  gravity  tends  to  stretch  the 
disturbance  wavelengths  as  the  axial  position  of  the  wave  increases. 

The  authors  gratefully  acknowledge  the  support  of  this  work  by  the  Air  Force  Office  of  Scientific  Research 
under  contract  number  F49620-92-J-0390. 


Nomenclature 

a  —  orifice  radius 

Cs  =  smoothing  constant 

D  =  doublet  contribution  matrix 

E(p)  =  complete  elliptic  integral  of  the  second  kind 

G  =  free  space  Greens  function 

Ii  —  IA  =  components  of  the  D  and  S  matrices 

K(p)  —  complete  elliptic  integral  of  the  first  kind 

hr  =  unit  normal  in  radial  direction 

hz  =  unit  normal  in  axial  direction 

Pn  =  Legendre  polynomial  of  order  n 

q  =  surface  velocity  in  the  normal  direction 

r  =  radial  coordinate 

S  =  source  contribution  matrix 

t  =  time 


z  =  axial  coordinate 

Bo  —  Bond  number,  Bo  =  pga2/<r 
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We  =  Weber  number,  We  —  pU2a/cr 
a  =  boundary  point  singularity  contribution 
f3  =  surface  slope 
T  =  domain  boundary 
k  —  surface  curvature 
p  =  density 
<r  =  surface  tension 
—  velocity  potential 
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Figure  7:  Typical  Droplet  Shedding  History,  We=1.45,  Bo=0.204  (equal  time  intervals). 
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Figure  8:  Droplet  Size  for  Dripping  Flow,  Bo=0.204. 
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Figure  11:  Jet  Profile  Near  Pinch  Location  Showing  Presence  of  Main  and  Satellite  Droplets,  We  =  100,  Bo 
=  0,  and  3%  Inflow  Velocity  Perturbation.  Equal  Time  Intervals  (dt=5.0),  with  the  First  Profile  at  t  =  87.00. 
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